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PREFACE 


This  report  was  prepared  by  the  D^artment  of  Civil  and  Enviromnental  Engineering,  West 
\Trginia  University,  P.O.  Box  6103,  Morgantown,  WV  26506-6103.  This  research  project  was 
funded  by  the  United  States  Air  Force  under  contract  F0863794M6010  with  the  325  Contracting 
Squadron,  Tyndall  Air  Force  Base,  FL  32403-5323.  The  Princ^al  hivestigator  was  Dr.  Donald 
D.  Gray.  The  point  of  contact  was  Dr.  Thomas  Stauffer,  OL-AL/EQ.  The  work  was  conducted 
between  March,  1994,  and  December,  1995. 

This  final  r^ort  describes  the  use  of  public  domain  software  to  create  numerical  models  of  the 
groundwater  flow  and  contaminant  tran^ort  in  the  Macrodn^ersion  E?q>eriment  2  conducted  at 
Columbus  AFB,  hfissis^>pi 

The  authors  thank  Dr.  Robert  EH,  Dr.  Mohammed  Gabr,  and  Mr.  Mandar  Sarangdhar  of  West 
\Trginia  University;  Dr.  Manfted  Koch  of  Ten]q)le  University;  and  Dr.  Chimmiao  Zheng  of  the 
University  of  Alabama. 
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EXECUTIVE  SUMMARY 


PabEc  domain  numerical  groundwater  flow  (MODFLOW)  and  contaminant  tran^ort 
(MT3D)  coniqputer  programs  were  used  to  simulate  groundwater  flow  and  contaminant  migration 
of  a  nonreactive  tracer  and  four  organic  coiiq>ounds  throu^  a  site  at  Columbus  Air  Force  Base  in 
Columbus,  Mississi^i  TEe  MADE-2  fMacrodispersian  Experiment)  experiment  commenced 
with  tile  pulse  injectum  of  the  dissolved  contaminants  into  the  saturated  zone  of  an  alluvial  aquifer 
through  a  series  of  five  injection  weDs.  The  concentrations  of  the  contaminants,  hydraufic  head, 
and  net  recharge  were  monitored  for  15  months.  Previous  measurements  of  hydrauEc 
conductivity,  poroaty,  and  other  parameters  were  also  used  in  the  modeling  process. 

The  finite>difference  code,  MODFLOW,  solved  the  governing  groundwater  equation  fin 
hydraufic  head  and  sewage  velocity  at  discrete  nodes.  Kiigmg  was  enqrloyed  to  obtain  inirial  and 
boundary  conditions  by  extnqrolating  measured  heads. 

The  head  and  flow  velocities  were  used  as  input  to  the  tran^ort  code.  MT3D  solved  the 
contaminaiit  tran^ort  equation  separate^  fiom  the  flow  equation  since  buoyancy  effects  were 
assumed  to  be  negligible.  Initial  snnnlations  of  the  tritmm  phune  were  unsatisfectory.  The 
simulated  plume  (fid  not  extend  fir  enough  fiom  the  injection  ate  to  match  the  observations.  An 
ad  hoc  assumption  of  horizontal  anisotropy  was  ^rpfied  to  the  conductivity  fidd  in  order  to 
increase  the  longitudinal  velocities  gnongh  to  pudi  the  plunie  downgraifient.  This  produced 
realistic  smmlated  plumes  for  trithim  and  tire  four  organic  ghemitsala 

Altiiough  the  shnnlated  plumes  were  realistic,  the  horizontal!^  anisotropic  conductivity 
field  used  to  achieve  these  pre(fictians  was  not.  The  assumed  piinc^al  axes  were  not  aligned  with 
the  directions  suggested  by  pump  tests  or  geologic  history,  hi  addition,  the  d^ree  of  anisotropy 
was  fir  in  excess  of  any  previous^  observed.  Fhialfy,  the  magmtnH#  of  the  effective  conductivity 
was  fir  above  \viiat  had  been  measured  at  the  site. 


V 


The  use  of  public  domain  sofiwaie  to  model  contaminant  tran^oit  at  heterogeneous  sites 
is  technical^  practical;  but  in  the  absence  of  exhaustive  field  measurements,  little  confidence  can 
be  placed  in  the  predictions. 
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SECTION  I 
INTRODUCTION 


A.  OBJECTIVE 

Niunerous  Air  Force  sites  have  suffered  groundwater  contamination  due  to  hydrocarbon 
^£05.  The  design  of  optimal  remediation  programs  requires  reliable  numerical  modeling  of 
contaminant  tranq>ort,  and  this  modeling  should  be  carried  out  using  the  most  widely  available 
hardware  and  software.  The  objective  of  this  stu^  was  to  illustrate  state-of-the-art  modding  of 
dissolved  contaminant  tran^ort  in  a  heterogeneous  suifidal  aquifer  using  public  domain  software 
on  a  microconDputer.  This  was  accon:plished  by  appfying  these  techniques  to  the  data  obtained  in 
an  actual  field-scale  experiment 

B.  BACKC»OUND 

Realizing  the  dangers  assodated  with  polhited  groundwater,  the  Air  Force  has  initiated 
several  programs  dedicated  to  understanding  the  physical  processes  involved  with  the  ftte  of 
contammants  in  the  snbsurfiice.  Focusing  mainfy  on  organic  solvents,  sudi  as  cleaming  solutions 
(e.g.  trichloroethylene)  and  jet  fuel  constituents  (e.g.  n^fathalae,  toluene,  xylene,  benzene,  etc.), 
these  studies  are  aimed  at  develophig  a  predictive  cq>abi]ity  for  the  tranq)ort  of  tirese  harmful 
chemicals.  By  estimating  the  amount  and  location  of  the  contaminant  at  any  given  time, 
remediation  techniques  can  be  hxptemented  in  a  more  cost-effective  maimer. 

In  tile  present  stiufy,  a  numerical  modd  was  appUsd  to  the  MADE  (Macrodispersion 
Experiment)  she,  located  at  Cohmdms  Air  Force  Base  (CAFB),  hfississqipi  Figure  1  show^  the 
location  of  the  site.  The  MADE  site  was  the  location  of  two  natural-gradient,  large-scale  tracer 
experiments  mtoided  to  quantify  tiie  effects  of  mamodi^ersivity  in  an  extremdy  heterogeneous 
aqui&r.  Macrodi^ersion  arises  ftom  spatial  variability  in  the  hydraulic  properties  of  natural 
aquifers  (Bo^  et  aL,  1992),  and  its  understanding  requires  controlled  e?q>erimentation  in  which  a 
detailed  anafyss  is  performed  on  a  nonreactive  tracer  plume.  Describing  the  motion  of  the  plume 
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requires  a  three-dimensicHial  saiqplmg  network  to  characterize  the  dq>th  and  areal  exteat  of 
contamination.  Extensive  hydraulic  conductivity  measurements  are  also  needed  to  accurate^ 
estimate  the  ^adal  variability  of  the  conductivity.  Together,  these  data  allow  the  diversion 
coefficients  to  be  estimated,  allowing  a  model  to  fit  the  observed  phenomena. 

The  MADE  site  is  the  most  heterogeneous  ate  used  to  date  for  a  natural-gradient 
experiment.  Table  1,  taken  firom  Reh&ldt  (1992X  shows  the  variance  in  the  natural  logarithm  of 
the  conductivity  (ln(K))  to  be  4.5,  translating  to  a  range  of  over  three  orders  of  magnitude.  This 
value  was  calculated  from  2187  separate  bordiole  flowmeter  measurements  in  49  diflerent 
profiles  by  Rehfeldt  et  aL  (1992). 


TABLE  1.  CONDUCUVriY  AT  MACRODISPERSION  STIES. 


K  =  hydraulic  conductivity  (cm/s) 

Lh  ==  horizontal  correktion  scale  (meters) 
Lv = vertical  correlation  scale  (meters) 


Location 

Variance  (hi(iQ) 

u 

L, 

MADE 

4.5 

12.8 

1.6 

Borden 

0.29 

2.8 

0.12 

Cape  Cod 

0.26 

5.1 

0.26 

Twin  Lakes 

0.031 

3.0 

0.91 
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Figme  1.  Location  of  MADE  Site.  Source:  Boggs  et  aL  (1993). 


Two  e?q)eiimeaQts  took  place  at  the  MADE  site,  denoted  as  MADE-1,  and  MADE-2.  The 
first,  MADE-1,  tracked  a  conservative  tracer,  bromide,  for  ^roximately  20  months.  The 
e^eiiment  started  with  a  pnlse  injection  of  a  known  quantity  of  solution,  and  concentrations  were 
monitored  periodical^  by  use  of  multilevel  saiq>lers.  The  resuhs  of  this  experiment  can  be  found 
in  a  series  of  journal  articles  (Boggs  et  aL,  1992;  Adams  and  Gelhar,  1992;  Rehfeldt,  Boggs,  and 
Gdhar,  1992;  Boggs  and  Adams,  1992). 

The  second  experiment,  MADE-2,  ^onsored  by  the  Air  Force  and  the  Electric  Power 
Research  histitute  (EPRI),  is  considered  in  this  research.  MADE-2  studied  the  transport  of  a 
conservative  tracer,  tritiated  water,  as  well  as  the  effects  of  sorption  and  biodegradation  on  four 
nonconsovative  dissolved  organic  compounds  (benzene,  naphthalene^  p-^ene,  and  o- 
dichlorobenzene)  in  the  saturated  zone.  The  dissolved  organic  chemicals  are  of  types  fi^und  in  jet 
fiids  and  cleaning  solvents.  MADE-2  started  in  lime,  1990,  wifii  the  pulse  injection  of  file 
solution  into  five  screened  wells,  spaced  1  meter  q>ait  The  infection  of  9.7  m^  of  solution  lasted 
48.5  hours.  For  IS  mcmfiis,  the  concentrations  of  eadi  contaminant  were  monitored  by  the 
Tomessee  Valley  Authority  (TVA)  using  multilevel  sanplers  di^ersed  somewhat  regularly  over 
the  domain.  The  mnltilevel  sampler  networic  provided  a  detailed  array  of  contaminant  plume  data. 
Five  three-dimensional  sanopling  events,  referred  to  as  “siuq>shots”  were  performed  at  intervals  of 
approximately  100  days.  The  features  of  the  observed  plumes  were  diaracterized  by  file  first 
three  ^atial  moments:  the  zeroth  (total  massX  first  (centroidX  and  second  (concentration  variance 
about  the  centroid)  moments. 

In  addition  to  file  extensive  concentration  data,  several  aquifer  parameters  were  measured 
to  determine  the  variability  of  file  hydraulic  conductivity  field.  First,  the  borehole  flowmeter 
measured  the  horizontal  hydraulic  conductivity  in  77  separate  boreholes  at  IS  an  (6  inch)  dqjth 
intervals.  Ifydranlic  conductivity  was  oflen  seen  to  range  over  fimr  orders  of  magnitude  in  the 
same  bordiole,  making  the  MADE  site  the  most  heterogeneous  aquifer  in  which  a  natural- 
gradient  tracer  experiment  has  been  conducted,  including  the  Cq>e  Cod  (LeBlanc  et  aL,  1991)  and 
Borden  sites  (Sudiclty,  1986).  Over  2500  measurements  allowed  the  use  of  geostatistical 
techniques  to  estimate  conductivity  values  over  the  entire  domain. 
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A  soil  grain  analyss  was  also  used  to  measure  the  hydraulic  conductivity.  The  empirical 
method  of  Seiler  (Boggs  et  al,  1990)  was  used  to  calculate  the  conductivity  from  grain-size 
distributions.  However,  onfy  214  soil  sauries  from  30  separate  coreholes  were  taken.  The 
variance  of  the  natural  logarithm  of  conductivity  measurements  calculated  to  be  3.1, 

outside  of  the  confidence  limits  of  the  bordiole  flowmeter  measurements,  3.4  to  5.6.  This  high 
degree  of  uncertainty  in  the  variance  estimates  for  the  borehole  flowmeter  is  unavoidable  for 
practical  problems  involving  extremely  heterogeneous  aquifers  (Rdifeldt  et  aL,  1992).  Other 
methods  used  to  measure  conductivity  at  the  MADE  site  were  the  slug,  peimeameter,  and  double¬ 
packer  tests,  having  variances  of  1.8, 5.5,  and  0.47,  respectivefy. 

HydrauHc  head  measurements  were  takoi  to  moxutm  the  rise  and  £dl  of  the  water  table.  A 
total  of  17  monthfy  surveys  were  taken  in  an  array  of  49  angle-  and  multistaged  piezometers. 
Continuous  hydrogr^hs  were  measured  in  8  staged  piezometers. 

C.  SCOPE 

Ihe  above  data,  supplied  by  the  Tennessee  Valley  Authority  and  tiie  United  States  Air 
Force,  were  used  to  spfify  and  evaluate  a  numerical  model  of  the  MADE-2  phnnes.  During  the 
Summer  Faculty  Research  Program  (SFRP)  of  1992,  Gray  selected  MODFLOW  (McDonald  and 
Haibaugh^  1988)  as  the  code  to  amnlate  the  thre^dhnensional  groundwater  flow  (Gray,  1992). 
Prehminaiy  steady  state  solutions  were  obtained  at  this  time. 

Gray  (1993)  continned  his  research  during  the  summer  of  1993  in  the  SFRP  at  Tyndall  Air 
Force  Base.  He  modeled  the  trithnn  phune  niong  the  mixed  Fnlerian-T  jgnmgim  finh-tvHiflferenfA 
program  MT3D,  ^^hich  solves  the  contaminant  tran^ort  equation  n-qng  the  MODFLOW- 
predicted  flow  fidd.  Two  transiient  flow  Adds  were  simulated,  one  assuming  a  unifi>nu 
conductivity  field,  and  the  otiier  a  heterogeneous  conductivity  field  based  on  kriging  the  borehole 
flowmeter  data.  The  uniform  conductivity  flow  modd  ran  approximate^  three  timas  fei^er  to 
simulate  a  given  time  period.  Ifowever,  the  MT3D  transport  simulations  were  unsuccessfiil 
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because  munerical  tnctahiKtieg  or  grossfy  unrealistic  predictions  ended  every  siinulation 
prematurely. 

Dining  the  gimmer  of  1994  at  Tyndall  Air  Force  Base,  Gray  and  Rucker  (1994)  were  able 
to  formulate  a  more  accurate  depiction  of  the  conductivity  field  by  using  a  new  geostatistical 
program,  Geo-EAS.  Again,  MODFLOW  and  MT3D  were  used  to  solve  the  flow  and 
concentration  fields,  re^ectively.  Their  successes  in  modeling  the  tran^ort  of  tritium  were 
limited,  as  die  code  ran  excessively  long.  Numerical  instability  again  ended  every  tran^ort 
simulation  before  the  derired  conq>letion  time. 

The  authors  continued  woric  on  the  MADE-2  modeling  project  under  the  present  contract. 
However,  it  was  not  until  the  ^ring  of  1995  that  the  tran^ort  of  tritium  was  modeled  for  the 
entire  468-day  experiment.  Excessive  run  times  and  inaccurate  results  did  not  allow  many 
rinmlations  to  be  corrpleted.  Tbe  tran^ort  code  took  iq>  to  3  weeks  to  finidi. 

Rucker  was  invited  for  a  second  summer  during  1995  to  continue  research  at  Tyndall  Air 
Force  Base.  A  modified  veraon  of  MT3D,  acquired  fixun  Dr.  Manfred  Koch  of  Tenq>le 
University,  reduced  run  times  from  weeks  to  hours.  Aidditional  modeling  investigated  the 
possibility  of  horizontal  anisotropy  in  the  conductivity.  Eatfy  successfiil  shnnlations  showed  the 
longitudinal  migration  of  tritium  to  readi  only  75  meters  downgradient  fix>m  the  injection  source. 
However,  tire  field  observations  tiiowed  the  plume  to  spread  at  least  225  m.  By  inqtlementing 
horizontal  anisotropy,  the  velocities  were  increased,  thus  increasing  the  phune  movement 
Continuing  tire  work  at  WVU,  the  rtithim  plume  was  finally  simulated  accurately,  and  other 
runs  were  conducted  to  smmlatethe  dissolved  organic  contaminants. 

In  addition  to  the  works  of  GnQ^  (1992,  1993)  and  Gray  and  Rucker  (1994),  there  have 
been  other  models  of  the  MADE-2  plumes.  Dr.  Manfred  Koch  conducted  independent  modding 
of  MADE-2  while  at  Tyndall  AFB.  Kodi  (1994)  was  able  to  model  the  contammants  through  the 
entire  468-day  e?q>etimait,  but  could  not  matdi  the  field  observations.  Dr.  C.  Zheng  and  Dr.  J. 
Jiao  of  the  University  of  Alabama,  and  C.  J.  NeviUe  of  S.  S.  Papadopulos  and  Associates  (Zheng 
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et  aL,  1994)  appfied  a  steady-state  model  to  the  MADE-2  e^erimoxt.  Their  efforts  proved  that 
the  simulated  phone  was  more  sensitive  to  the  hydraulic  conductivity  field  than  to  the  di^ersivity 
used  Each  study  greatly  contiibuled  to  the  present  research 

D.  APPROACH 

A  model  is  an  approximate  representation  of  an  observed  phenomenon.  A 
mathematical  model  consists  of  the  governing  equations  together  with  appropriate  boundary  and 
niitial  conditions.  The  mathematical  equations  may  be  solved  analytically  or  numerically,  with 
numerical  solutions  usually  allowing  for  more  complex  boundaries  and  heterogeneous  soil 
properties  (Anderson  and  Woessner,  1992).  This  study  concerns  numerical  modeling. 

The  partial  differential  equations  which  govern  groimdwater  flow  and  contaminant 
tran^ort  can  be  discretized  by  the  finite-difference  method,  the  finite-element  mi^od,  or  some 
other  schemes.  The  finite-difference  method  is  the  most  easify  understood  and  was  used  in  thfg 
study.  The  domain  is  giidded  into  regulaify  shaped  blocks  or  ceDs.  Within  each  cell  are  nodes 
^^ere  the  unknowns  are  to  be  calculated.  The  discretized  equations  are  solved  nging  numerical 
methods  implemented  in  conq)uter  programs  known  as  codes.  The  codes  used  in  this  research 
were  the  groundwater  flow  code  MODFLOW  and  the  contaminant  transport  code  MT3D.  Both 
are  widefy  used  public  domain  programs  \rixich  can  be  executed  on  current  personal  conqmters. 

This  report  provides  a  detailed  e?q>lanation  of  a  nmneiical  tnodel  being  appfied  to  simulate 
the  MADE-2  field  data.  The  r^ort  moves  ff om  the  site  descrqrtion  to  data  anafysis,  mnHftKng 
and  ends  with  the  interpretation  of  results  and  suggestions  fin  furtiier  research. 
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SECTION  n 
SITE  DESCRIPTION 

A.  HYDROGEOLOGIC  MEASUREMENTS 

Figure  1  shows  that  the  MADE-2  test  site  is  located  approximately  6  km  east  of  the 
Tombigbee  River  and  2.5  km  south  of  the  Buttahatchee  River  at  Columbus  Air  Force  Base  in 
Lowndes  Coimty,  Mississippi-  The  site  lies  within  the  Cohimbus  aquifer  and  is  situated  above  the 
100-year  flood  plain.  Averaging  8  to  10  meters  in  depth,  the  Cohimbus  aquifer  is  con^osed  of 
alluvial  deposits,  primaifly  inteifliigeiing,  discontinuous,  lozenge-sh^ed  lenses  of  pooify  to  well 
sorted  gravefy  sand  and  sandy  gravel  with  small  amounts  of  clay  and  silt  (Young,  1994).  The 
lenses  are  typica%  on  the  order  of  8  meters  in  the  horizontal  direction,  and  1  meter  in  the  vertical 
Subsurfece  investigations  show  the  soil  to  be  mainfy  unconsofidated  and  cohesionless,  however 
occasional  small  consolidated  zones  were  encountered  during  eiqiloratoiy  drilling  (Boggs  et  aL, 
1992).  Beneath  the  afluvial  aquifer  lies  the  Eutaw  fbimation,  an  aquitard  consisting  of  clays,  flne 
grained  sands,  and  silts  of  marine  origin,  estimated  to  be  75  meters  thick. 

Data  gathered  at  the  MADE  site  mduded  hydraulic  conductivity,  hydraulic  head,  aquifer 
porosity,  and  bulk  density.  These  extensive  measurfements  allowed  anafysis  of  the  ^atial 
distribution  of  hydraulic  conductivity  using  geostatistical  techniques.  Moreover,  the  data  have 
been  used  to  infer  the  geologic  history  of  tire  area. 

L  Hydraulic  Conductivity 

Ifydraulic  conductivity  was  measured  using  several  direct  and  indirect  metitods.  Of  the 
direct  methods,  indudhig  bortiiole  flowmeter,  permeameter,  sing,  and  double-padter  tests,  the 
bordiole  flowmetCT  tests  provided  the  largest  data  set,  with  over  2500  measurements  in  77 
different  locations.  Bordiole  flowmeter  measurements  were  originally  made  fer  the  MADE-1 
eT^eriment,  but  21  new  weDs  have  been  tested  since  that  time.  The  borehole  flowmeter  measures 
the  discharge  during  punq)ing  at  discrete  locations  along  tire  screened  length  of  the  wdL  Ibis 
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provides  the  inflow  to  the  well  as  a  function  of  elevation.  Since  the  inflow  is  assumed  to  be 
horizontal  and  radial,  the  disdiarge  firom  each  layer  is  proportional  to  the  product  of  the  layer 
hydraulic  conductivity  and  the  layer  thickness.  Using  the  layer  flow  rate  and  drawdown,  the 
Cooper- Jacob  well  equation  is  employed  incrementally  to  evaluate  the  hydraulic  conductivity  of 
each  layer.  The  key  assunq>tions  underlying  this  procedure  are  ( 1)  the  aquifer  is  layered  and  each 
layer  is  homogeneous  and  of  uniform  thickness,  (2)  the  storage  coefficient  of  each  layer  is  hneariy 
related  to  the  layer  transmissivity,  and  (3)  the  well  losses  attributed  to  each  layer  can  be  estimated 
(Boggs  et  aL,  1990).  The  heterogeneity  of  the  MADE  site  Hmits  ffie  validity  of  these 
assunq>tions,  e^edaDy  assunq>tion  1. 

Flow  measurements  were  made  with  an  impeller-type  flowmeter  which  was  lowered  down 
the  bordhole.  The  rotation  of  the  inq>eller  caused  by  vertical  flow  in  the  well  was  read  by  optical 
sensors  and  converted  to  a  voltage  whidi  was  directly  proportional  to  the  rate  of  rotation 
(Rehfeldt  et  aL,  1989).  Ihe  voltage  passed  through  a  coaxial  cable  to  surfece  electronics. 
Calibration  of  the  mstmment  was  necessary  before  and  afler  each  flowmeter  measurement  to 
convert  the  recorded  voltages  to  flow  rates.  Ihe  lower  thre^old  of  flow  measurement  for  the 
inq>eller  flowmeter  is  approximate  O.OOS  Us,  coiTeq>onding  to  a  hydraulic  conductivity  of  about 
10“*  cm/s  (Rdifeldt  et  aL,  1989). 

The  locations  of  the  bordiole  flowmeter  measurements  can  be  seen  in  Figure  2.  The 
origin  of  the  plot  lies  at  the  center  of  the  injection  wells.  Each  conductive  is  labeled  with  a 
‘K’  and  a  number  corresponding  to  foe  order  in  which  it  was  installed.  The  weDs  were 
constructed  of  S.l  cm  diameter  fln^  PVC  dotted  p^e,  and  were  soeaied  over  the  entire 
saturated  lengffi  of  the  aquifer,  excqrt  for  gaps  where  sections  were  joined  (Rehfeldt  et  aL,  1992). 
Within  each  weD,  flowmeter  measurements  were  made  at  {^proximate  15.24  cm  (6  inch)  intervals. 
For  a  more  detailed  ejqilanation  of  the  borehole  flowmeter  method  or  the  derivation  of  the  model 
equations  used  to  calculate  the  hydrauHc  conductivity,  see  Rdifeldt  et  aL  (1989). 
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Figure  3  shows  the  dqith-averaged  hydraulic  conductivity  variation  and  a  profile  along  a 
vertical  slice.  The  injection  site  lies  within  an  area  of  relatively  low  hydraulic  conductivity  (10'^ 
cm/s),  and  the  conductivity  increases  one  to  two  orders  of  magnitude  to  the  nortL  At  the 
northern  extreme  of  the  domain  the  mean  conductivity  decreases  to  values  noticed  in  the  near 
field.  The  portion  of  the  aquifer  above  57  meters  MSL  shows  a  higher  conductivity  variation  than 
the  lower  portion.  A  region  of  hi^  conductivity  from  southwest  to  northeast  throu^  the 
midsection  of  the  aquifer  can  also  be  seen. 


The  areas  of  hi^  permeability  within  the  aquifer  have  been  attributed  by  Young 
(1994;199S)  to  a  former  meandering  river  charmel,  developed  during  the  end  of  the  Pleistocene 
period.  Among  the  evidence  cited  by  Young  (1995)  fin  this  hypothecs  is  an  aerial  photograph  of 
the  test  site  (Flgtne  4)  r^ttch  shows  a  dififexence  in  v^etation  crossing  the  site  at  about  30**  east 
of  north.  The  charmd  is  said  to  fie  above  58  meters  MSL  and  is  about  70  meters  wide, 
corie^onding  to  the  zone  of  higher  conductivity.  The  majrn  argument  used  by  Rdb&Idt  et  aL 
(1992)  to  iuterpr^  this  band  as  nothing  more  fiian  a  regional^  continuous  zone  of  relative^  high 
mean  conductivity  is  that  the  observed  groundwater  flow  does  not  seem  to  fiillow  the  axis  of  the 
channel,  but  moves  perpendicolariy  to  the  northwest.  These  considerations  will  be  seen  to  have  a 
significant  bearing  on  the  mndelmg  process. 
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Figure  3.  Vertically  Averaged  I^draolic  Coadncdvity  and  Candacdvity  Profile  Along  Section 
A-A.  Source:  Boggs  et  aL  (1993). 
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Figure  4.  Aerial  Fbotogi^h  of  MADE  Site. 


Laboratoiy  penneam^er  tests  of  recompacted  soil  sainples  raised  questions  of  whether  the 
measured  hydraulic  conductivity  was  representative  of  the  in*situ  conductivity  due  to  the 
cohesionless  nature  of  the  highly  sandy  soils.  Tests  on  mhiimally  disturbed,  trimmed  samples  from 
lined  ^ht  core  barrel  sanq>lers  were  similariy  su^ect.  The  solution  was  die  design  of  a  constant 
head  penneabilrty  test  which  could  be  conducted  without  removing  the  soil  from  the  sampling 
tube.  Ei^ty-eight  measurements  were  taken  from  nine  core  samples.  The  permeameter  and 
borehole  flowmeter  measurements  differed  by  iq>  to  three  orders  of  magnitude  in  some  areas. 
This  discrepancy  was  attributed  in  part  to  the  fret  that  the  permeameter  measures  vertical 
conductivity  wfrereas  the  borehole  flowmeter  measures  horizontal  values.  A  further  reason  is  the 
difference  in  effective  measurement  volumes  between  the  two  instrumoits.  An  additional 
limitation  of  the  permeameter  data  is  that  the  nine  core  samples  came  from  the  .same  area  and 
cannot  characterize  the  entire  site. 

Slug  tests,  which  consist  of  instantaneous^  adding  a  known  volume  of  water  to  a  weh  and 
recording  the  decline  in  pressure  head  over  time,  were  conducted  in  22  partially  penetrating 
piezometers  in  order  to  measure  hydranfic  conductivity.  Inconclusive  results  were  obtained,  as 
the  undersong  assunqitions  of  the  mathematical  models  were  violated.  The  extreme  heterogeneity 
of  file  aquifer  conqiromised  the  tests  agnifiicantfy,  as  rqiorted  by  Boggs  et  aL  (1990). 

The  double-packer  test  consists  of  tiie  iogection  of  water  into  an  isolated  interval  of  the 
weQ  under  constant  pressure  or  constant  flow  rate.  The  isolation  is  achieved  by  the  inflation  of 
packers  located  above  and  below  the  iiyection  interval  The  results  of  the  packer  tests  gave 
horizontal  hydranfic  conductivities  fiur  discrete  layers  that  were  consistent^  higher  than  those 
measured  witii  the  borehole  flowmeter  method.  Boggs  et  aL  (1990)  concluded  that  the  di^aiity 
is  due  to  artificial  vertical  movement  of  the  hqected  water  within  the  disturbed  ammhs  of  soil  just 
outside  of  the  well  caang.  Aiiq>le  evidence  points  to  the  disturbance  of  sediments  a<yacent  to  the 
well,  Mhich  may  have  increased  the  hydranfic  conductivity  conqiared  to  undistuibed  sediments. 
The  flowmeter  method  is  less  susceptible  to  annulus  effects,  because  punping  of  the  fiifiy 
screened  well  produces  predominantly  horizontal  flow  towards  the  well  (Boggs  et  aL,  1990). 
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The  indirect  calculations  of  hydraulic  conductivity  included  grain-size  analysis,  direct- 
current  resistivity,  and  sur&ce  geophysical  surveys.  The  follo\ving  discussion  ^  be  hmited  to 
the  soil  grain-size  analysis  as  the  others  are  beyond  the  scope  of  this  report  Boggs  ^  aL  (1990) 
discussed  the  estimation  of  the  hydraulic  conductivity  from  enq)irical  formulas  using 
measurements  of  grain-size  distributions.  A  total  of  214  measurements  were  made  on  soil 
sanq>les  located  q>oradica%  around  the  ate.  The  geometric  mean  of  the  hydraulic  conductivity 
was  much  higher  than  fotmd  by  the  direct  methods. 

An  additional  study,  conducted  by  Stauffer  and  Manoranjan  (1994),  proved  by  vertical 
kiiging  and  segmented  trend  surfeces,  the  extent  to  vdiich  grain-size  analy^  can  be  usefiiL  Core 
sauries  coUected  over  the  entire  saturated  length  of  the  aquifer  from  17  irregularly  ^aced 
locations  were  dried  and  sieved  to  determine  particle  size  distributions.  The  enqnrical  formula 
used  in  the  preceding  experiment  by  Boggs  et  al  (1990)  was  again  employed  to  determine 
hydraulic  conductivity.  However,  the  data  was  not  analyzed  in  two  or  three  spatial  dimeasions  to 
determine  variability  in  the  data.  On  contrary,  the  data  was  kriged  in  onfy  one  dimension, 
vertically,  to  obtain  statistical  insi^  into  the  trends  in  conductivity.  Urey  concluded  that  grain- 
size  anafysis  can  provide  understanding  of  the  general  trend  of  hydraulic  conductivity  patterns 
vsduch  are  comparable  to  those  found  with  the  bordiole  flowmeter. 

2.  Aquifer  Tests 

Two  large-scale  traditional  aquifer  tests,  titled  ATI  and  AT2,  were  conducted  at  the  site 
to  estimate  the  average  hydraulic  conductivity  and  ^edfib  yield.  The  first,  ATI,  started  hr 
March,  1985,  and  lasted  ea^  days.  For  three  days  drawdown  was  measured  in  twelve  primary 
observation  wells  by  pressure  transdncers  and  a  data  logging  system.  The  observation  weDs  were 
positioned  in  three  radial  lines  leading  fimn  the  punqnng  weO,  PWl,  whidi  lay  outside  the  domain 
of  the  MADE-2  plumes,  approximatdy  70  meters  southeast  of  tire  hyection  wells.  Thirteen 
additional  wells  within  a  lOO-nwter  radius  of  PWl  were  also  monitored  twice  daily  dmmg  the  test 
in  order  to  delineate  the  extent  of  the  drawdown  cone  (Boggs  et  aL,  1990). 
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The  period  of  punoping  was  followed  by  5  days  of  recovery  measuiements.  The  resultant 
plots  of  time  versus  drawdown  were  evaluated  using  the  Neumann  type-curve  method  for 
anisotropic  confined  aquifers.  From  the  plots,  the  mean  transmisstvity  and  ^edfic  yield  were 
calculated  as  1.8  cmVs  and  0.04,  re^ectivdy.  Assuming  an  average  saturated  thickness  of  8.2  m, 
the  mean  hydraulic  conductivity  was  0.002  cm/s.  In  addition,  the  drawdown  curves  from  each 
observation  well  were  used  to  construct  drawdown  contours  whidi  were  distorted  aliases, 
suggestive  of  horizontal  heterogeneity  or  anisotropy. 

AT2  was  conducted  within  the  MADE-2  experimental  area,  with  pumping  well  PW2 
located  65  meters  north  of  the  injection  ate.  Observation  wells  were  aligned  along  three  radials 
extending  fiom  PW2.  A  total  of  32  staged  observation  weDs  (20  primary  and  12  secondary)  were 
monitored  for  the  16  day  test  (8  days  of  pumping  and  8  days  of  recovery).  The  different  stages  of 
the  weDs  allowed  anafysis  of  the  vertical  anisotropy  of  conductivity.  Hie  drawdown  data  were 
plotted  on  log-log  grq>hs  to  obtain  transmissivity  and  ^edfic  yield  values  using  the  Neumann 
type-curve  method.  Tlie  average  transmissivity  value  of  the  test  was  calculated  by  aiitiuneticalfy 
averaging  all  values  obtained  fiom  the  pumping  well-observation  well  pairs.  The  transmissivity 
value  of  20.1  cm^/s  was  divided  by  the  average  saturated  thickness  of  9.8  meters  to  obtain  a 
hydraulic  conductivity  of  0.02  cm/s.  These  values  are  an  order  of  magnitude  larger  than  in  ATI 
because  of  the  diflferent  location,  which  is  mote  rqnresentative  of  the  MADE-2  plnme  region. 
The  mean  specific  yield  in  AT2  was  d^ermined  to  be  0.1. 

Contour  plots  of  the  drawdown  were  eOptical  as  seen  with  ATI.  It  is  ddiated  whether 
heterogeneity  or  horizontal  anisotropy  is  the  cause.  A  value  of  2.6  was  conputed  for  the 
horizontal  anisotropy  with  the  Tnavmnrm  pnncpal  axis  oriented  about  35**  west  of  north.  A 
vertical  anisotropy  (  Ky/EQ,  )  was  also  calculated  fiom  the  mean  of  the  staged  observation  weDs; 
the  value  was  0.18. 

Table  2  summarizes  the  average  hydranfic  conductivity  values  obtained  fiom  each  test 
Calculations  for  the  average  conductivity  in  the  aqni&r  tests  only  included  the  primary  wdDs. 
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TABLE  2.  SUMMARY  OF  HYDRAULIC  COM)UCTIVnY  MEASUREMENTS. 


Method 

No.  of  Measurements 

Mean  Conductivity  cm/s 

Borehole  Flowmeter 

2187 

5.52  X  10’’ 

Permeameter 

88 

6.12x10'^ 

Slug  Test 

22 

1.65  X  10‘* 

DoublO'Packer 

37 

4.20  X  10'^ 

Grain-Size  Analysis 

214 

4.30x10’^ 

ATI 

11 

2.00  X  10'^ 

AT2 

15 

2.00  X  10'^ 

3.  Hydraulic  Head  Monitoiing 

Mooitoimg  of  the  hydraulic  (piezometcic)  head  was  conducted  before  and  doting  the 
tracer  e7q)erinient  using  angle  and  multistage  piezometers.  Figure  5  shows  the  bcadons  of  the 
piezometers.  The  weDs  are  designated  with  a  ‘P’  and  a  number  corresponding  to  the  order  in 
vsbidi  they  were  placed.  A  sofBx  of  ‘A’  ,  ‘B’,  or  both  follows  the  piezometer  number, 
designating  the  level  of  screening  of  the  weR  An  *A’  r^resents  soeoung  of  the  iq>per  level;  a 
‘B’  denotes  die  lower  level;  the  combination  indicates  multilevel  screening.  The  actual  screen 
levels  varied;  but  lower  stage  was  screened  at  an  average  elevation  of  56.3  meters,  and  the  iqtper 
at  an  average  elevation  of  61.1  meters.  Forty-nine  piezometers  were  monitored  manually  using  an 
electric  probe  on  a  monthfy  basis.  A  total  of  17  mondify  surveys  were  taVai^  starting  June  19, 
1990,  and  ending  on  S^tember  11,  1991.  Moreover,  eig^  pahs  of  staged  piezometers  were 
equ^ed  with  contnmous  water  level  recorders  (piezometers  P8AB,  P27AB,  P44AB,  P53AB, 
P54AB,  P55AB,  P60AB,  P61AB).  Measurements  taken  &om  the  condnuous^  monitored 
piezometers  matched  the  nMimal  survey  heads  clos^. 
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The  groundwater  table  fluctuated  seasonally  over  a  range  of  2-3  meters,  resulting  in  an 
approximate  30  percent  variation  in  saturated  thickness  of  flie  aquifer  (Boggs  et  aL,  1993).  The 
average  horizontal  hydraulic  gradient  of  0.0017,  ^erienced  seasonal  periodicity  corre^onding 
to  water  table  fluctuations  (Stauffer  et  aL,  1994).  Figure  6  shows  hydraulic  heads  from  Jime  19, 
1990,  for  both  iq>per  and  lower  screened  piezomet^  Heads  were  kriged  using  <$0-EAS  and 
plotted  using  the  contouring  and  suifrice  mapping  program  SURFER  5.0  for  Windows.  The 
hydraulic  heads  d^  towards  the  plan  northwest.  In  the  southern  end  the  contours  are  closely 
packed,  reflecting  the  lower  values  of  transmissivity  in  the  southern  zone  conq)ared  to  the  mid  and 
frr  field  wfrere  the  head  contours  are  more  ^read  out 
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Hydrogrq>lis  were  produced  widi  the  data  obtained  from  the  continuously  monitored 
piezometers  showing  head  versus  days  since  injection.  The  hydraulic  heads  were  averaged  over 
approximate  monthly  intervals  for  coizg)arison  with  computer  simulations  described  later  in  the 
paper.  Figure  7  shows  a  plot  of  two  averaged  hydrographs  from  piezometers  located 
approximately  50  meters  ^art  Well  P53  is  located  in  an  area  of  low  conductivity,  5  meters 
iq)gradient  from  the  injection  point.  Well  PS4,  located  20  meters  downgradient  from  the  injection 
site,  is  in  an  area  of  hi^er  conductivity.  There  are  about  two  orders  of  magnitude  difference 
between  the  two  locations.  Ihe  cross-correlation,  computed  for  the  two  data  sets  by  the 
FORTRAN  code  CROSS,  written  by  Dr.  Wilson  of  West  \%ginia  University,  is  presented  in 
Figure  8.  The  figure  shows  a  lag  of  0  (zero)  days,  indicating  a  simultaneous  rise  in  head. 
However,  when  the  cross-correlation  is  computed  for  two  piezometers  separated  by  a  distance  of 
approximately  280  meters  (piezometer  P53  and  piezometer  P61,  located  264  meters 
downgradient  fix>m  injection),  the  lag  is  -1  days.  The  results  inq>]y  that  the  rise  in  head  in 
downgradient  well  P61  h^ens  approximately  one  day  earlier  than  in  P53,  decreasing  the 
hydraulic  gradient  in  times  of  hi^er  recharge  to  the  aquifor.  A  possible  justification  for  this 
phenomenon  is  that  P61  lies  in  an  area  of  lower  conductivity  compared  to  PS3.  Recharge  thus 
produces  a  small  mounding  efEect  in  P61,  before  the  groundwater  level  reaches  an  equifibrium 
However,  the  lag  has  a  minimal  effect  on  foe  flow,  and  it  appears  that  foe  water  table  tends  to  rise 
almost  simultaneous  over  foe  entire  domain. 
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Cross  Correlation  of  Piezometers  P53  and  P54 


Lag  (days) 

Cross  Correlation  of  Piezometers  P53  and  P61 


Lag  (days) 


Figure  8.  Cioss-Coirdadon  of  Piezometers  P53  and  P54  (Upper)  and  Piezometers  P53  and  P61 

(Lower). 
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4.  Other  Aquifer  Measurements 


Aqui&r  porosity  and  bulk  density  were  calculated  from  84  mininiany  disturbed  soil  core 
samples  collected  from  four  separate  core  holes  at  the  tracer  test  site  (Boggs  et  aL,  1992).  The 
average  porosity  of  the  aqui&r  was  calculated  to  be  0.3 1  with  a  standard  deviation  of  0.08.  The 
trend  in  the  porosity  showed  higher  values  as  dq)th  increased.  The  lower  elevations,  below  about 
56  meters  MSL,  exhibited  values  of  0.4  and  greater,  while  the  ipper  levels  decreased  to  as  low  as 
0.2.  This  is  consistent  with  the  soil  types  found  in  the  repecdve  devations:  sand>filled  gravel 
with  open  or  clay>filled  por^  in  frte  ipper  layers  and  a  hi^er  sand  content,  averaging  70%,  in  the 
lower  levels. 

The  dry  bulk  density  was  calculated  from  the  same  core  samples  as  the  porosity 
measurements.  The  volumes  of  the  core  samples  were  measured.  The  samples  were  then  oven 
dried  and  sieved  to  find  particle  density.  The  bulk  denaty  was  then  estimated  from  the  particle 
density  and  file  measured  volume.  The  average  value  for  the  bulk  density  was  1.77  g/cm^  wifii  a 
standard  deviation  of  0.18  g/cm^ 

B.  HYDROLCKjY 

Daily  precpitation  and  temperature  data  were  collected  at  the  CAFB  weather  station, 
located  less  than  2  km  away,  to  quantify  the  effects  of  recharge  to  the  Columbus  aquifer.  Daify 
pan  evaporation  data  were  collected  at  Mississip^  State  Univasity,  approximately  35  km  distant, 
and  supplied  by  the  State  Climatologist,  Dr.  C.  L.  Wax.  Based  on  the  recommendation  of  Dr. 
Wax,  a  pan  coefficient  of  0.8  was  used  to  estimate  the  evapotranspirafion  (Gray  and  Rucker; 
1994).  The  net  reduuge  was  then  calculated  by  subtracting  an  estimated  evpotranspiration  value 
from  daily  precpitation.  Missing  evaporation  data  were  estimated  from  the  daily  maxhnnm 
temperatures  using  the  empirical  equatian  of  Pote  and  Wax. 
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A  gr^h  of  the  net  rechaige  versos  time  (days  since  injection)  is  presented  in  Figure  9. 
Negative  values  in  the  graph  depict  days  that  had  higher  evapotran^iration  than  precipitation. 
These  are  more  common  during  the  summer  months. 
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SECTION  m 
•ERACER  MONITORING 


The  tracers  for  MADE'2  were  tiitiated  water  and  four  dissolved  organic  confounds 
(benzene,  nq)hthalene,  p-?Qiene,  and  o-dkhlorobenzene).  Tritium  was  used  as  a  passive  tracer 
since  it  does  not  undergo  chemical  reactions  or  sorption  and  has  a  half  life  much  longer  than  the 
period  of  the  e?q)eriment.  Eie  four  organic  corr^ounds  are  common  constituents  of  fiiels  and 
solvents.  A  small  fiacdon  of  the  p-xylene  was  labeled  with  radioactive  to  study  the  biological 
or  chemical  transformation  of  p-xylene  should  it  occur  duriog  the  e?q>etinient.  Table  3  lists  the 
initial  concentrations  and  total  injected  mass  for  each  tracer. 

The  expoiment  began  with  the  injection  of  9.7  m^  of  solution  into  five  weDs,  ^aced  1 
mrter  apart  in  a  linear  array.  The  wells  were  screened  from  57.5  to  58.1  meters  MSL.  The 
injection  started  on  June  26,  1990,  and  lasted  48.5  hours.  A  constant  iryection  rate  of  3.3  Urma 
was  maintained,  raising  the  hydraulic  head  in  the  injection  wehs  by  0.45  meters.  The  contaminant 
solution  was  prqpared  and  stored  on  site  using  ambient  groundwater  from  a  well  located  iq)stream 
from  the  injection  point.  No  non-aqueous  phase  cmitaminants  were  injected. 


TABLES.  INITIAL  CONCENTRATIONS  AND  INJECTED  MASS. 


Tracer 

Mean  Concentration 

Mass  Injected 

Tritium 

55610  pCi/mL 

0.5387  a 

■Hi3iS!SS1iH 

2770pCi/mL 

0.0268  a 

benzene 

68.1  mg/L 

659.7  g 

p-xylene 

41.4  mg/L 

402.0  g 

naplithalene 

7.23  mg/L 

70.0  g 

o-didilorobenzene 

32.8  mg/L 

317.7  g 

The  monitoring  of  the  tracers  was  accoiqpfished  by  withdrawing  groundwater  samples 
through  multilevel  saiiq)lers  (MLS)  and  poative  di^lacement  (BarCad)  sanqilers.  Figure  10, 
ad^ted  from  Boggs  et  al  (1993),  shows  the  placement  of  the  MLS  and  BarCad  samplers.  A 
total  of  328  MLS  were  placed  in  the  study  area  in  a  complex  ^atial  pattern.  Each  MLS 
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incorporated  20  to  30  san^Hng  points  ^aced  0.38  mieters  apart  in  the  vertical  directi<m.  A  three 
dimensional  rq>resentation  of  the  phune  was  in&rred  fiom  the  MLS  data. 

A  total  of  five  MLS  sanqiling  episodes,  or  “sn^shots”  were  conducted  at  approximate 
100  day  intervals,  starting  27  days  after  injection.  Snapshots  1-4  surveyed  the  entire  phime;  but 
Snapshot  5,  starting  440  days  after  injection,  was  designed  onfy  to  investigate  the  bounds  of  the 
organic  plumes  and  did  not  encompass  the  more  extoistve  tritium  or  plumes. 

Additional  san^ling  of  the  contaminaats  was  accoiiq)fi^ed  using  BarCad  positive 
di^lacement  saiq>lers,  placed  along  two  “fencehnes”  oriented  normal  to  the  flow.  Ihe  ftncelines 
were  two  parallel  rows,  approximate^  6  and  16  miners  downgradient  firom  the  hyection  ate. 
Sampling  took  place  every  two  weeks  in  the  earfy  stages  of  the  e^eriment,  eventually  moving  to 
3-month  intervals  in  the  latter  portion.  Table  4  lists  tiie  tracer  sampling  dates  of  both  sn^shot 
and  finceline  data  with,  the  total  immber  of  weDs  sampled  at  each  time. 
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(Ill)  A 
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Flgme  10.  Saiiq>fiiig  WeUNetwodc.  Source:  Boggs,  et  aL  (1993). 
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TABLE  4.  SUMMARY  OF 'mACER  SAMPLING. 


Sampling  Event 

(F=fencelmef 

S=snapshot) 

Date 

Time  After  Injection 
(days) 

No.  of  Wells  Sampled 

FOl 

Jul  9-11, 1990 

13 

26 

S21 

Jul  23-27,1990 

27 

99 

F02 

Aug  13-17, 1990 

48 

31 

F03 

Sep  17-19, 1990 

83 

53 

F04 

Oct  5-17,1990 

111 

39 

S22 

Nov  5-8, 1990 

132 

111 

F05 

Dec  3-4,1990 

160 

29 

F06 

Jan  8-9, 1991 

195 

25 

S23 

Fd)  5-7,1991 

224 

190 

F07 

281 

42 

S24 

May  21-23,  1991 

328* 

205 

S25 

Sep  9-11, 1991 

440 

79 

*  May  21,  1991,  is  actualfy  329  days  since  the  stait  of  ngection.  Boggs  et  aL  (1993),  Machityre 
et  aL  (1993),  and  Stauffer  et  aL  (1994)  all  refer  to  May  21  as  day  328.  That  convention  is 
maintained  in  this  report 


Analyses  of  tiitmm  and  were  conducted  at  the  Water  Resources  Research  Center  at 
hfisas^>pi  State  University.  The  field  saiiq>les  were  measured  with  a  liquid  scintillation  counter 
in  dual-isotope  mode.  The  background  concentrations  of  tritium  and  at  tiie  site  established 
the  sensitivity  of  the  measurements  and  woe  finmd  to  be  2  and  3  pCi/mL,  respective^. 


Analytical  measurements  were  performed  on  the  organic  tracers  by  the  TVA 
Envirmnnental  Chenastiy  Laboratory.  The  extracted  organic  tracers  were  anatyzed  by  gas 
chromatogr^hy  (GC)  using  a  flame  ionization  detectmr  (FID)  tystem.  The  sensitivity  for  the 
GC/FID  method  was  4  )i%/L  for  naphthalene,  p-x^ene,  and  o>DCB,  and  50  |rg/L  for  benzoie. 
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SECTION  IV 
MADE-2  DATABASE 


Die  database  for  MADE-2,  conqiiled  by  the  Tennessee  Valley  Authority  (TVA),  is  stored 
on  three  3.5-inch  disks.  DISKA  contains  the  concentration  data  for  benzene,  n^hthalene,  p- 
xylene,  o-dichlorobenzene,  tritimn,  and  DISKB  stores  the  piezometiic  head  measuronents, 
and  DISKC  stores  hydraulic  conductivity  data  from  the  borehole  flowmeter  tests.  Additional 
information,  including  format  and  organization,  can  be  found  in  the  READ.ME  files  located  on 
each  disk. 

The  tracer  concentration  data  are  separated  into  two  directories  :  \SNAPSHOT  and 
\FENCE.  Under  \SNAPSHOT  are  located  the  five  thiee-dimiensional  plume  snapshots.  Each 
sn^shot  is  contained  in  individual  £Qes  mariced  SNAP#.DAT,  \^ere  #  is  the  sn^shot  number. 
Table  5  is  an  exanqile  of  the  format  used  fi^  the  snqi^ot  data. 


TABLES.  EXAMPLE  SNAPSHOT  DATA. 


BflQZBie 

p~Xyimt 

oDCB 

Trithm 

Cttboft-14 

X 

Y 

Z 

DaysSmoe 

Sas^lelD 

(^) 

(ugflL) 

(pCML) 

(pCML) 

(m) 

(m) 

(m) 

Date 

Ixjectkn 

S22D00606 

660 

63 

550 

720 

1697 

45.7 

-13 

67 

5662 

11/5/90 

132 

S22D00608 

1800 

260 

1800 

1700 

3464 

140.9 

H 

67 

57.12 

11/5/90 

132 

S22D00610 

4500 

480 

4200 

3100 

7024 

264 

-13 

67 

11/5/90 

132 

S22D0a614 

2000 

200 

1800 

1500 

3277 

129.7 

m 

67 

5665 

11/5/90 

132 

S22D00616 

3800 

400 

4716 

205.4 

H 

67 

jggj 

11/5/90 

132 

S22D00S1Z 

2900 

290 

2600 

3715 

158L1 

m 

67 

59.66 

11/5/90 

132 

$22000620 

jjjgjU 

320 

2900 

2200 

4222 

169.4 

-13 

67 

60.17 

11/5/90 

132 

$22000622 

1300 

67 

590 

450 

2714 

9&8 

m 

67 

61.08 

11/5/90 

132 

S22D008Q3 

52 

<  4.0 

9 

HI 

3 

33 

43 

42 

55.84 

11/5/90 

132 

S22D00805 

88 

<  4.0 

12 

8 

57 

5.1 

43 

42 

561 

11/5/90 

132 

S22D00807 

63 

<  4.0 

27 

36 

143 

73 

43 

42 

56.85 

11/5/90 

132 

S22D008Q9 

81 

5 

61 

64 

41 

53 

43 

42 

5736 

11/5/90 

132 
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Under  the  directoiy  \  DISKA  \  FENCE  are  located  the  tracer  data  associated  with  the 
fenceHne  sanq>lers.  The  format  is  similar  to  the  snapshot  data  and  a  s^arate  file  is  designated  for 
each  data  set. 

The  piezometer  data  on  DISKS  is  divided  into  two  subdirectories,  \MONTHLY,  and 
VRECORDER.  The  “monthly”  surveys  are  found  under  \MONTHLY,  and  a  separate  file  is 
provided  fi>r  survey.  H^meter  labels  ending  in  “A”  are  screened  in  the  i^per  level,  and  labels 
ending  in  “B”  designate  lower  level  screens.  Table  6  gives  an  fflcasqple  of  the  piezometiic  data. 


TABLE  6.  EXAMPLE  SURVEY  PIEZOMETRIC  HEAD  DATA 


wdl 

iill 

date 

dapsed  time 
(days) 

wato*  level  devation 
(meters) 

-86.1 

2.56 

1/8/91 

196 

63.82 

P-40 

-11.41 

83.28 

1/8/91 

196 

63.47 

P-41 

1.62 

56.17 

1/8/91 

196 

63.51 

P-45 

-30.84 

10.38 

1/8/91 

196 

63.77 

P-52 

-85.85 

190.51 

1/8/91 

196 

63.05 

P-8A 

95.13 

123.75 

1/8/91 

196 

r  63.49 

P-8B 

95.13 

123.74 

1/8/91 

196 

63.41 

P-lOA 

103.2 

5.73 

1/8/91 

196 

65.43 

P-lOB 

101.69 

6.15 

1/8/91 

196 

63.73 

Hie  subdirectory  VRECORDER  of  DISKS  contains  the  continuousfy  monitored 
piezometer  data  for  16  weDs.  An  example  of  the  data  is  provided  in  Table  7. 

TABLE?.  EXAMPLE  CONTINUOUS  PIEZOMETRIC  HEAD  DATA 


wdl 

z 

(meters) 

y 

(meters) 

dapsedtime 

(days) 

year 

Julian  day 

water  levd  devation 
(metars) 

P-53a 

-4.9 

-10.3 

-7 

90 

170 

62.44 

P.53a 

-4.9 

-10.3 

-6 

90 

171 

62.42 

P-53a 

-4.9 

-10.3 

-5 

90 

172 

62.4 

P-53a 

-4.9 

-10.3 

-4 

90 

173 

62.39 

P-53a 

-4.9 

-10.3 

-3 

90 

174 

62.36 

P-53a 

-4.9 

-10.3 

-2 

90 

175 

62.33 

P-53a 

-4.9 

r  -10.3 

-1 

90 

176 

62.31 
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The  hydraulic  conductivity  data  for  67  of  the  weOs  tested  with  bordiole  flowmeters  are  on 
DISKC.  Each  well  has  a  separate  flle  designated  by  the  well  name.  The  remaining  10  profiles 
were  measured  afler  the  conchision  of  MADE-2,  and  the  data  are  available  from  TVA  on  a 
separate  disk.  The  authors  have  incoiporated  these  new  data  on  their  own  copies  of  DISKC. 
Table  8  is  an  exanople  of  a  conductivity  data  file. 


TABLE  8.  EXAMPLE  BOREHOLE  FLOWMETER  CONDUCTTVITY  DATA. 


FLOWMETER  WELL  K-14  HYDRAULIC  CONDUCTTVITY  -  DEPTH  PLOT  DATA 
WELL  COORDINATES:  X=  114.42  Y  =  204.20 
[A]  DEPTH  BELOW  GRADE  (M) 

|B]  DEPTH  BELOW  GRADE  (FT) 

[C]  ELEVATTON  ABOVE  MEAN  SEA  LEVEL  (M) 

[D]  ELEVATION  ABOVE  MEAN  SEA  LEVEL  (FT) 

[E]  HYDRAUUC  CONDUCTTVITY  (CM/SEC) 


[F]  HEAD  IN  THE  AQUIFER  (FT  ABOVE  SEA  LEVEL) 


lAl 

m 

[C] 

m 

m 

4.654 

15.27 

60.67 

199.05 

8.30E-02 

201.4 

4.807 

15.77 

60.518 

198.55 

8.30E-02 

201.4 

4.807 

15.77 

60.518 

198.55 

3.23E-02 

201.4 

4.959 

16.27 

60.366 

198.05 

3.23E-02 

201.4 

4.959 

16.27 

60.366 

198.05 

3.32E-02 

201.4 

5.111 

16.77 

60.213 

197.55 

3.32E-02 

201.4 

The  MADE>2  coordinate  system  (XY)  has  its  origin  at  the  center  of  the  injection  site  with 
the  longitudinal  axis  (Y-axis)  aligned  along  file  expected  mean  trqectoiy  of  file  plume,  12°  west 
of  norfiL  Since  borehole  flowmeter  data  were  earned  over  fiom  MADE-1,  many  of  the  wells  are 
stin  located  in  the  MADE-1  reference  coordinate  system  (XY^).  The  MADE-2  origin  is  located 
at  coordinates  (X’  =  85.2  meters,  Y’  =  188.4  meters)  wifli  the  Y-axis  rotated  25.68® 
counterclocJcwise  from  the  Y'-axis  of  MADE-1.  The  transformation  can  be  performed  by 
applying  the  following  formulas  on  wells  KOI  through  K59: 
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X’  =  Xcose-Ysine  +  h 
Y’=Ycos0-Xsiiie  +  k 

X  =  CX’-h)cos0  +  (Y’-k)sme 
Y  =  (Y’-k)cose  +  (X’-h)sm0 

where  X’  and  Y’  desdgnate  the  MADE-1  coordinates,  X  and  Y  are  the  MADE-2  coordinates,  h 
85.2  meters,  k  =  188.4  meters,  and  0  =  -25.68°  (Boggs  et  aL,  1993). 
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SECTION  V 
FLOW  MODELING 


A.  CODE  SELECTION 

A  numerical  code  had  to  be  chosen  to  model  MADE-2.  At  the  direction  of  die  ^(msor, 
only  public  domain  codes  available  at  Tyndall  Air  Force  Base  at  the  time  of  selection  (June  1992) 
were  considered.  Table  9  is  an  inventory  of  die  candidate  groundwater  codes  (Gray  ,1992). 

All  codes  Hsted  in  Table  9  solve  the  groundwater  flow  equation  and  the  advective- 
di^ersve  tran^ort  equation,  except  for  MODFLOW,  \^ch  requires  an  additional  code  (MT3D) 
to  solve  the  tran^ort  problem.  Qnfy  MODFLOW,  HST3D,  and  SWICHA  can  solve  the 
equations  in  all  three  spatial  dimensions.  Though  limited  to  Qn&>dimensional  problems, 
SAMFTID  can  predict  the  motion  of  up  to  three  immiscible  phases;  but  the  others  are  single 
phase  codes.  The  entry  ‘Unsat?’  refers  to  the  ability  of  the  code  to  solve  the  flow  ecpiation  in  the 
unsaturated  zone.  This  is  a  difBcult  task,  since  the  hydiauhc  conductivity  of  the  vadose  zone  is  a 
fimction  of  the  degree  of  saturation.  SUTRA  and  SAMFTID  have  this  c^abilitjr;  the  others  are 
onfy  valid  in  the  saturated  zone.  Pre-  and  posQnocessors  (denoted  as  ‘Pre?’  and  ‘Post?’)  are 
available  fiir  some  of  the  codes  to  allow  the  user  to  more  easify  manipulate  data  for  preparation  of 
hput  files  and  diph^r  of  ouput  files.  While  not  strict^  necessary,  these  programs  are  extreme^ 
usefuL 
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TABLE  9.  AVAB-ABLE  (mOUNDWATER  CODES  IN  JUNE  1992. 


Modd,  Version 

Date 

Flow 

X- 

port 

Dim 

Unsat? 

Pre? 

Post? 

Method 

MODFLOW,  4.2 

11/91 

yes 

no  * 

3 

no 

yes 

yes 

FD 

HST3D,  1.5 

2/92 

yes 

yes 

3 

no 

no 

no 

FD 

SWICHA,  5.05 

2/91 

yes 

yes 

3 

no 

no 

no 

FE 

SUTRA,  0690- 

2D 

6/90 

yes 

yes 

2 

partial 

yes 

yes 

FD/FE 

MOC,  3.0 

11/89 

yes 

yes 

2 

no 

yes 

no 

FD/MOC 

Random  Walk 

81 

yes 

yes 

2 

no 

no 

no 

FD/RW 

SAMFTID,  1.0 

9/90 

yes 

yes 

1 

yes 

? 

? 

FD/FE 

1 

H)  =  Finite-DifEeraDice,  FE  =  Fiiiite>Ekmeot,  MOC  =  Method  of  Characteristics 
RW  =  Random  walk 

*  Conq>anion  tran^ort  code  (MT3D)  now  available. 

In  view  of  the  extreme  heterogoidty  of  the  aquifer  and  die  nature  of  the  observed  plume, 
it  was  obvious  that  a  three-diniensional  solver  would  be  needed  to  produce  realistic  predictions 
(Gray,  1992).  From  the  codes  listed  in  Triile  9,  MODFLOW  was  diosen  fiir  its  ease  of  use, 
excellent  documentation,  and  wide  acceptance. 

MODFLOW  (modular  thre^dimensional  fimte-difiEerence  flow  modd)  was  written  by 
McDonald  and  Harbaugh  (1988)  of  the  U.  S.  Geological  Survey.  Original^,  MODFLOW  was 
coded  in  FORTRAN  66,  but  was  upgraded  to  FORTRAN  77.  MODELX>W  has  a  modular 
structure,  \dierein  similar  program  functions  are  groiq)ed  together,  and  ^edfic  computational 
and  hydrologic  options  are  independent  of  other  options.  Such  structming  allows  the  addition 
and  subtraction  of  new  modules  as  the  need  arises  without  the  diai^on  of  the  rest  of  the  code. 
The  major  options  can  simulate  the  effects  of  wdls,  recharge,  rivers,  drains,  evapotranspiration, 
streams,  and  general  head  boundaries.  Ihe  solution  methods,  which  solve  the  matrix  equations 
established  by  the  MODFLOW,  are  found  in  the  strong^  inplicit  module  (SIP),  the  slice- 
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successive  over  relaxation  module  (SOR),  and  the  preconditioned  conjugate  gradient  module 
(PCG2),  developed  by  IM  (1991). 

The  input  structure  of  the  program  uses  separate  batch  files.  Format  instructions  stated 
Avithin  the  batdi  file  dictate  the  format  for  the  input  without  modification  to  the  program.  The 
type  of  output  may  also  be  selected  to  fit  a  particular  need.  An  ou^ut  control  module  allows  the 
results  to  be  saved  on  didc  or  printed  to  the  screen. 

MODFLOW  solves  the  groundwater  flow  equation  (1)  for  the  hydrauhc  head  using  a 
finite-difference  approximation.  Equation  (1)  assumes  time-dependent,  constant  density 
groundwater  flow  in  an  anisotropic,  heterogeneous,  saturated  medium. 


\^here:  x,y,  z  are  the  principal  coordinates. 

K»(,K^,Kg2  are  the  piincpal  hydraulic  conductivity  values. 

h  is  the  hydraulic  (piezomebic)  head. 

W  is  the  net  volumetric  inflow  or  outflow  per  volume  of  aquifer  (sources  /  sinks). 

S,  is  the  specific  storage. 

Bxploying  either  head  or  flow  boundary  conditions  and  initial  (starting)  heads  with 
Equation  (1)  constitutes  a  mathematical  rq>resentation  of  a  groundwater  flow  system  (McDonald 
and  Harbaugii,  1988).  The  solution  ^es  head  values  as  a  fimction  of  pace  and  time.  Specific 
discharges  are  estimated  by  difEerenong  tile  heads.  The  assunption  of  constant  density  inpficit  in 
Equation  (1)  means  that  buoyancy  forces  are  neglected  and  decoiples  tiie  tranport  equation  firom 
the  flow  equation.  This  permits  one  to  solve  the  flow  problem  before  considenng  the  tranport 
problem  The  validity  ofthis  assumption  will  be  discussed  later. 
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B.  DISCRETIZAITON 


To  iiiq>leniieiit  a  numerical  model,  a  proper  grid  must  be  chosen.  The  MADE-2  site  is  an 
area  of  ^roximat^  300  meters  x  200  meters  with  about  2  meters  of  relief  The  dimioisions  of 
the  conqmtational  grid  were  reduced  to  330  meters  by  105  meters,  tising  a  uniform  grid  pacing 
of  5  meter  x  5  meter  cells.  Grid  convergence  tests  will  be  discussed  later.  MODFLOW  can 
simulate  flow  uang  a  vaiiabfy'sized  grid,  however  the  uniformity  allows  smq>lificatiQn  in  the 
extraction  of  results  and  aquifer  param^ers  (hydraulic  head,  permeability,  ^edfic  yield,  etc.). 
The  saturated  zone,  vriiose  thickness  varied  feom  9  to  11  meters,  was  discretized  into  9  layers,  as 
described  below.  66  rows  and  21  columns,  the  number  of  cells  totaled  12  474. 

In  terms  of  MADE>2  coordinates,  the  domam  stretched  feom  -52.5  meters  to  52.5  meters 
in  the  X  direction  and  -27.5  meters  to  302.5  meters  in  the  Y  direction,  with  the  axes  parallel  to 
the  MADE-2  coordinate  ^em.  Ahhou^  the  Y  axis  runs  12"*  west  of  true  north,  the  positive  Y 
direction  is  defined  as  plan  north.  The  directions  mentioned  in  later  sections  of  this  rq>ott  are 
plan  directions  unless  otherwise  ^ecified.  The  origin  of  the  domain  was  the  site  of  the  central 
injection  well  AH  five  injection  wells  were  contained  in  one  cell  (Row  61,  Column  1 1).  Figure 
11  shows  the  grid  used  fbr  the  flow  model,  and  later  fin  the  tran^ort  model 
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Figure  11.  Gnd  Used  in  MADE-2  Modeling.  Note  Vertical  Exaggeration. 
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MODFLOW  uses  a  block-centered  tedmiqae  to  evaluate  the  conductances  between 
ai^acent  nodes.  The  original  block-centered  flow  package  (BCFl)  can  shuulate  both  confined  and 
unconfined  aquifer  systems.  It  allows  cells  to  dry  out  as  the  water  table  but  cannot  allow 
rewetting  of  cells  as  the  water  table  rises.  This  Kmitgrion  forced  an  earlier  discretization  of 
MADE-2  to  \ise  a  thick  top  layer  to  insure  that  the  top  layer  never  went  dry  (Gray,  1993).  In  that 
grid,  the  bottom  eight  layers  were  1  meter  thick  and  the  top,  with  a  base  at  59.0  meters  MSL,  was 
as  much  as  6  meters  thick.  A  later  version  of  the  block-centered  flow  package  (BCF2),  developed 
by  McDonald  et  al  (1991),  permits  rewetting  as  the  water  table  rises  above  the  bottom  of  a  dry 
cdL  This  is  an  inq)ortant  feature  for  simulating  the  MADE-2  e3q>eriment,  as  the  water  table  was 
seal  to  fluctuate  between  2-3  meters. 

The  discretization  enqiloyed  in  this  work  set  the  base  of  the  top  layer  at  63.0  meters  MSL, 
so  that  the  saturated  thickness  of  the  top  layer  never  exceeded  2. 1  meters.  The  next  seven  layers 
were  each  1  meter  thidc;  and  the  bottom  layer  varied  flx>m  56.0  meters  MSL  at  its  top,  to  the 
impermeable  bottom  formed  by  the  aquitard  bdow.  The  thickness  of  the  lower  layer,  averaging 
3.31  mders,  ranged  from  0.49  meters  at  the  thinnest  point  to  6.1  meters  at  the  thickest  This  is  a 
large  rang^  but  the  bottom  layer  is  the  least  hiq^ortant  to  the  overall  flow.  In  terms  of 
MODFLOW  classification,  layer  1  is  unconfined,  layers  2  through  7  are  fiilly  converiible,  and 
layers  8  and  9  are  confined  (Gray  and  Rucker,  1994). 

Tenqioral  discretization  divided  the  468-day  oqierimeDt  into  stress  periods  and  time  stqis. 
Stress  periods  are  defined  in  regard  to  MODFLOW  as  time  intervals  in  wdiich  all  external  stresses 
(sources  and  boundary  condidoins)  are  constant  (McDonald  and  Ibtban^  1988).  In  turn  the 
stress  periods  are  divided  into  time  stq>s.  The  stress  periods  for  the  MADE-2  experiment  were 
centered  on  the  17  piezometiic  surveys  and  the  injection  period  firr  a  total  of  18  stress  periods. 
Eadi  stress  period  was  divided  into  2-day  time  steps.  Except  for  the  injection  period,  the  stress 
periods  are  roughfy  centered  on  tiie  head  survey  dates.  Table  10  lists  the  survey  dates  and  stress 
period  lengths  used  fin  the  MADE-2  shunlations.  Simulation  Day  numbers  are  counted  from  the 
start  of  the  shnulation  period  with  Sinmlation  Day  1  being  June  12,  1990.  The  injection  took 
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place  oa  Simulation  Days  15  and  16.  The  reader  should  be  aware  that  reports  presenting  field 
results  are  usualfy  givea  in  terms  of  days  since  the  start  of  injection. 


TABLE  10.  STRESS  PERIODS  FOR  MADE-2  SIMULATIONS. 


Stress 

Period 

Starting 

Date 

Starting 
Sim.  Day 
Number 

Stress  Period 
Length 
(Days) 

Head  Survey 
Date 

Head  Survey 
Sim.  Day 
Number 

1 

June  12, 1990 

1 

14 

June  19, 1990 

8 

2* 

June  26 

15 

2 

3 

June  28 

17 

36 

42 

4 

53 

28 

Aug.  13 

63 

5 

81 

32 

98 

6 

Oct  2 

113 

26 

Oct  15 

126 

7 

Oct  28 

139 

24 

Nov.  7 

149 

8 

Nov.  21 

163 

32 

Dec.  5 

177 

9 

Dec.  23 

195 

32 

Jan.  8, 1991 

211 

10 

Jan.  24, 1991 

227 

30 

Feb.  8 

242 

11 

F*.23 

257 

28 

March  8 

270 

12 

Marcii  23 

285 

30 

April  4 

297 

13 

315 

24 

333 

14 

339 

18 

343 

15 

Jime3 

357 

24 

June  13 

367 

16 

June  27 

381 

34 

393 

17 

415 

32 

■C!!Sfli 

434 

18 

Sept  1 

447 

22 

Sept  11 

457 

Last  day 

Sept  22 

468 

- 

m 

- 

*  hqecdon  period. 


C.  INmAL  AND  BOUNDARY  CONDITIONS 

The  piezometric  heads  firom  the  monthfy  surveys  were  used  to  establirii  the  initial  head  at 
each  node,  as  wdl  as  the  head  at  eadi  boundary  node  as  a  fimction  of  time  (Gray  and  Rocker, 
1994).  Because  piezometers  were  not  placed  at  every  nod^  the  data  were  kriged  to  infer  die 
needed  values.  Kriging  is  a  geostaristical  procedure  by  uhich  a  rdativety  qwall  number  of 
irregularfy  ^aced  data  are  used  to  estimate  values  at  a  large  number  of  discrete  points  using  a 
weighted  moving  average  interpolation  method.  Kriging  is  the  best  finear  unbiased  estimator  and 
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reproduces  the  measuremeots  exactfy.  Ihe  kriged  heads  from  the  first  survey  were  used  to 
egtahligh  the  heads,  with  the  kriged  results  of  the  lata  surveys  being  used  as  boundary 
conditions  in  MODFLOW’s  General  Head  Boundary  module.  Earlier  modeling  by  Gray  (1993) 
assumed  the  hydraulic  head  did  not  vary  as  a  fimction  of  depth  and  enoployed  Surfer  Version  4  to 
krige  heads  in  two  ^atial  dimensions.  These  were  used  as  boundary  and  initial  conditions  for 
every  layer.  However  Gray  and  Rucker  (1994)  recognized  that  the  heads  were  a  fimction  of 
depth.  They  also  used  a  new  geostatisdcal  program,  Geo*£AS,  which  allowed  more  flexibility  in 
the  krigmg  process. 

Geo-EAS  Version  1.2. 1,  developed  by  Enghmd  and  Sparks  (1991)  for  the  EPA,  is  a  menu 
driven  geostatisdcal  program  whidb  performs  two-dimensional  kriging.  It  allows  the  user  to  plot 
the  variogram  and  change  tire  variogram  to  fit  the  data.  Curroitfy,  titree  variogram  types  can  be 
used  in  Geo-EAS:  linear,  q>heiical,  and  exponential  Eadb  incorporates  a  nugget  and  sill  to 
propedy  reflect  the  data’s  variation.  Surfor  Version  4,  used  by  Gray  (1993X  does  not  give  the 
flexibility  that  Geo-EAS  allows.  The  variogram  onty  uses  a  linear  model  whh  the  efifects  of  the 
nugget  and  sQl  being  overiooked.  However,  it  does  produce  better  contour  plots  of  the  results. 
Later  versions  of  Surfor  have  removed  many  of  the  dtortcormogs  mentioned  ^ove. 

Hie  commercial  ^readdieet  Quattro  Pro  was  used  to  evaluate  the  distribution  of  the 
piezometer  screm  mi^oints.  Hie  average  elevations  of  the  iq>per  and  lower  piezometer  screens 
were  found  to  be  60.5  and  56.0  meters  MSL,  re^ectivety.  Based  on  the  average  screen 
elevation,  the  piezometers  were  divided  into  two  groiq>&  For  exanple,  the  file  GW010891 
(groundwater  head  survey  of  January  8,  1991)  was  divided  to  create  UP018091.DAT  and 
L0010891J>AT.  Hie  program  MADETCKX  refinmatted  the  files  to  Geo-EAS  input 
requirements  and  «^Knwn  piezometers  v^duch  were  for  outside  the  computational  domain  or 
were  not  screened  close  to  the  average  elevations.  A  total  of  15  piezometers,  vdiose  outpoints 
ranged  from  59.76  to  61.22  meters  MSL,  were  used  in  the  upper  set  of  files.  The  lower  set  of  23 
piezometers  had  mu^oints  v\d)ich  ranged  from  55.51  to  56.71  meters  MSL.  Hie  horizontal 
distribution  of  the  upper  and  lower  piezometers  seen  in  Figure  12  shows  ^arse  coverage  towards 
the  north  end  of  the  giil 
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Geo-EAS  was  used  to  sqpaiatefy  kiige  the  upper  and  lower  piezometer  files  for  each 
survey  date.  Table  11  shows  the  parameters  used  in  fitting  models  to  the  variograms,  including 
model  type,  nugget,  ^  and  range.  With  a  &w  exceptions,  the  linear  model  gave  a  good  fit  to  the 
data.  As  anomalies  were  found  in  the  data,  piezometers  were  removed  to  allow  a  better  fit  For 
example,  the  file  UOP081991.DAT  had  an  odd  point,  with  the  piezometer  P43A  having  a  head  of 
58.78  meters  MSL.  Hiis  point  was  removed  from  the  variogram.  The  result;  a  better  looking 
variogram  with  the  curve  fitting  more  closed  to  the  data.  Figures  13  and  14  show  contour  m^s 
of  kiiged  monthly  survey  heads  for  June  19,  1990,  and  March  8, 1991.  The  contour  m^s  for  the 
remaining  ktiged  head  surveys  can  be  found  in  Appendix  A.  In  almost  every  survey,  the  heads  (% 
towards  the  northwest  The  contour  plots  were  cneated  by  Surfor  Version  5  for  Windows.  The 
program  (j£02SRF  reformatted  the  Geo-EAS  output  for  input  to  Surfer. 
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TABLE  11  (A).  KRIGING  PARAMETERS  FOR  THE  UPPER  LEVEL  HEAD  SURVEYS. 


TABLE  11  (B).  KRICHHG  PARAMETERS  FOR  THE  LOWER  LEVEL  HEAD  SURVEYS. 
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For  iiqtiit  to  MODFLOW  the  kiiged  heads  of  the  iq>per  piezometers  were  assigned  to 
Layers  1  tfaiou^  4.  The  kriged  heads  of  the  lower  piezometers  were  assigned  to  Layers  8  and  9. 
Linear  interpolation  was  used  to  assign  heads  to  Layers  5,  6,  and  7.  Program  BASMAKER  was 
written  as  a  preprocessor  for  MODFLOW  in  order  to  set  19  the  ‘basic’  input  data  inchiding  grid 
dimensions,  number  of  stress  periods,  hydrologic  packages  to  be  used,  and  mirial  head  at  each 
node.  Vertical  interpolation  of  the  initial  heads  was  performed  within  the  progranL 

The  boundary  conditions  for  the  model  were  established  in  the  General  Head  Boundary 
package.  This  package  served  to  q>ecify  heads  (a  Dirichlet  boundary  condition)  at  boundary 
nodes,  and  to  change  riiem  for  each  stress  period.  The  program  C^BMAKER  created  the  input 
rile  necessary  for  MODFLOW  execution  using  the  remaining  head  surveys.  The  program  extracts 
the  heads  onfy  along  the  boundary  and  uses  the  vertical  interpolation  scheme  described  above  to 
assign  boundary  heads  to  aU  the  layers. 

D.  MODELING  OF  AQUIFER  PARAMETERS 

Aquifor  parameters  such  as  horizontal  hydraulic  conductivity,  vertical  conductance,  and 
^edfic  yield  were  modeled  for  hqmt  to  the  Block-Centered  Flow  package.  Each  is  eqplained 
bdow. 


1.  Ifydraufic  Conductivity  and  Vertical  Leakance 

Horizontal  hydraulic  conductivity  was  measured  by  the  borehole  flowmeter  in  77  sq>arate 
profiles,  located  around  die  MADE-2  ate.  The  data  were  measured  over  successive  15-cm  layers 
within  die  saturated  zone  of  the  aquifer.  The  profiles  contained  g^s  caused  by  joints  in  the  well 
screens  ^Mdch  were  filled  in  using  values  feom  immediatdy  above  and  below  the  g^. 

In  order  to  combine  the  profiles  to  characterize  the  entire  site  it  was  inqiortant  to  note  that 
the  top  devadons  of  the  profiles  varied  from  57.62  meters  to  62.68  meters  MSL  dqieading  on-the 
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local  water  table  elevation  on  the  date  of  measurement.  One  consequence  was  that  the  15-cm 
layers  differed  from  one  profile  to  the  next.  KAVG94  was  written  to  average  the  borehole 
flowmeter  conductivities  over  the  MODFLOW  grid  layers.  TTie  program  extended  the  profiles  up 
to  64.0  meters  MSL  by  assuming  a  constant  conductivity  from  the  top  of  the  well  screen  to  the 
top  of  the  aquifer.  Tire  lowest  data  varied  from  5 1.88  meters  to  56.22  meters;  these  values  were 
extended  down  to  the  next  lower  integer  elevation-  The  15-cm  intervals  were  arithmetically 
averaged  over  each  MODFLOW  layer  to  generate  horizontal  conductivities  for  all  nine 
MODFLOW  layers  (Gray  and  Rucker,  1994).  This  process  assumes  horizontal  isotropy. 

Vertical  conductivity  betweai  the  layers  was  also  calculated  with  KAVG94,  assuming  the 
aquif^  material  was  locally  isotropic.  By  harmonically  avera^g  the  conductivity  between 
MODFLOW  nodes,  a  vertical  leakance  was  generated.  Vertical  leakance,  called  VCONT  in 
MODFLOW,  is  the  vertical  conductivity  divided  by  the  thickness  between  adjacent  nodes. 
MODFLOW  uses  VCONT  to  calculate  vertical  flow  between  successive  layers.  Because  of  the 
variable  thickness  of  the  lowest  layer  (9),  the  leakance  between  Layers  8  and  9  was  based  on  the 
interval  between  56.5  meters  and  55.5  meters,  except  for  three  profiles  (K-2,  K-26,  and  K-28) 
vriiich  ended  at  56.0  meters.  For  the  lowest  layer,  VCONT  is  inq)licitly  set  to  zero  because  the 
lower  boundary  of  the  domain  is  assumed  impermeable.  TTie  general  formula  for  VCONT  is 

VCONT^,.^  =  Av,  ^  Av>.. 

wbiere:  Av^  and  Avfc<.i  are  the  thicknesses  of  layers  k  and  k+1,  respective 
(Kz)ij^  and  are  the  vertical  conductivities 

A  major  disadvantage  of  VCONT,  as  seen  in  the  MADE-2  experiment,  happens  in  the  top  layers 
of  the  aquifer.  When  the  top  layer  is  unconfined  its  thickness  fluctuates  as  the  water  table 
The  average  water  table  depth  should  be  used  to  calculate  VCONT  between  the  water 
table  layer  and  the  one  below  it.  However,  VCONT  is  calculated  by  using  the  distance  between 
nodes,  since  the  water  table  elevation  is  not  known  a  priori  If  the  vertical  conductivity  is  not 
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homogeneous,  then  the  equivalent  conductivity  between  the  nodes  wfll  be  underestimated, 
because  the  iq)per  layer  is  assumed  thicker  than  it  really  is.  This  is  not  a  problem  for  the  top  most 
layer  (Layer  1)  as  it  is  always  homogeneous  (from  extrapolation  of  the  profiles  to  the  top  of  the 
aquifer).  But  the  problem  arises  in  lower  layers  as  the  water  table  reaches  them.  The  phenomena 
also  affects  horizontal  conductance,  winch  also  is  based  on  fully  saturated  conditions. 

The  output  file  for  KAVG94  (KAVG94.0UT)  contained  77  conductivity  profiles,  each 
averaged  over  the  same  nine  grid  layers.  The  well  number,  X  and  Y  well  coordinates,  grid  layer, 
horizontal  conductivity  (m/d),  transmissivity  (m^/d),  vertical  conductivity,  and  VCONT  were  all 

listed.  KAVTOGE  separated  the  output  into  layer  files  in  Geo-EAS  format  The  file  names  were 

K1AY#.DAT,  where  #  represented  the  MODELOW  layer  number.  The  layer  files  contained  only 
the  well  number,  location,  transmissivity  (conductivity  for  Layer  1)  and  VCONT.  In  order  to 
avoid  unphysical  negative  values  after  kriging,  and  to  req>ect  die  lognormal  distiibutimi  of  the 
conductivity,  the  data  were  log-transformed  by  the  program  KA2LOG  to  establish  ln(K)  values. 
The  data  were  presented  in  Geo-EAS  format  and  named  KLOG#.DAT.  Spherical  and 
e?q>onential  variograms  were  successfiiHy  fitted  to  the  log  transformed  data. 


The  data  for  each  layer  were  kriged  horizontalfy  to  obtain  natural  log  transmissivity  and 
natural  log  VCONT  values  at  every  grid  node.  The  data  were  then  transformed  back  to  unlogged 
values  by  the  program  DLOCHTLE.  The  files  were  converted  to  Surfer  format  and  contour  plots 
were  made  for  all  nine  layers.  Figure  15  shows  an  exanqile  of  the  kriged  transmissivity  and 
VCONT  for  Layer  3.  The  figure  shows  a  series  of  hi^  conductivity  regions  trending  from 
southwest  to  northeast  in  the  zone  from  about  Y  =  40  meters  to  200  meters.  These  may  reflect 
the  hypothesized  river  channeL 
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The  kiiged  conductivities  were  prepared  for  the  MODFLOW  input  file  of  BCF2  by  the 
program  BCF2MAKEL  This  preprocessor  is  project-specific,  Le.,  designed  for  the  ejqificit  use  of 
the  MADE-2  project.  The  program  also  added  the  layer  elevations,  horizontal  anisotropy  values, 
and  the  ^edfic  yield  (or  storage  coefficient  for  confined  aquifers)  to  the  input  file. 

2.  Horizontal  Anisotropy 

A  value  for  horizontal  anisotropy  can  be  introduced  in  the  MODFLOW  input  file.  Defined 
as  a  column  to  row  ratio,  the  value  multiplies  the  conductance  between  nodes  along  rows  to 
obtain  conductance  between  nodes  along  a  cohimn.  In  the  standard  version  of  MODFLOW  only 
one  value  can  be  entered  for  each  layer,  although  this  was  later  modified,  as  described  later.  As  a 
starting  point,  the  aquifer  was  assumed  to  be  horizontally  isotropic;  a  value  of  1.0  was  used  m  all 
nine  layers  for  the  initial  simulations. 

3.  Storage  Parameters 

In  a  steady-state  simulation,  the  ri^t  hand  side  (R.H.S.)  of  Equation  (1)  is  zero  (5h/a  = 
0)  and  qiedfic  yield  is  not  needed.  As  discussed,  the  experiment  underwent  tenqioral  changes  m 
the  hydrological  features,  requiring  a  transient  simulation.  Thus,  qiedfic  yield  values  were  needed 
to  properly  solve  the  groundwater  fl^ow  equation.  The  ^ecific  yield  was  diosen  based  on  the 
pumping  test,  AT2.  The  base  value  of  0.1  was  assigned  to  all  unconfined  layers  (including  the 
convertible  layers  as  they  dianged  firom  confined  to  unconfined).  No  measuremaits  were  made 
for  ^edfic  storage,  so  a  confined  storage  coefficient  base  value  of  0.0001  was  assumed,  based 
on  the  textbook  values  for  qiedfic  storage  in  sand  and  sandy  gravel  given  by  Anderson  and 
Woessner  (1992).  This  value  was  assigned  to  all  the  nodes  in  each  layer. 
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E.  HYDROLOGIC  STRESSES 


The  remaining  stresses  to  the  flow  system  include  recharge  and  injection  well  flow  rates. 
Since  the  site  was  covered  primaiify  by  weeds  and  brusii  and  no  sur&ce  water  was  observed,  the 
MODFLOW  sur&ce  water  packages  were  not  used. 

1.  Recharge 

The  daify  net  recharge  values  were  incorporated  into  the  flow  model  using  the  Recharge 
package.  Since  pan  evaporation,  used  to  account  for  evjq)otranq)iration,  was  subtracted  fi’om  the 
daily  prec^itation,  a  negative  value  of  recharge  indicated  a  net  loss  of  water.  Thus  the 
Evapotranq)iration  package  was  not  used.  The  daily  net  recharge  values  were  arithmetically 
averaged  over  each  stress  period  giving  the  values  seen  in  Table  12. 


TABLE  12.  AVERAGE  NET  RECHARC®. 


Stress  Period 

Stress  Poiod 

1 

-0.00313 

10 

0.00809 

2 

-0.00478 

11 

0.00114 

3 

-0.00148 

12 

0.00794 

4 

-0.00409 

13 

0.01022 

5 

-0.00286 

14 

0.00357 

6 

-0.00107 

15  1 

0.00046 

7 

0.00071  ^ 

16 

-0.00273 

8 

0.00942 

17 

-0.00159 

9 

0.00387 

18 

-0.00384 

2.  Well  Shnulation 

The  of  a  wcll  was  required  during  the  second  stress  period,  in  vriiidi  the 

injection  of  tracers  in  the  aqueous  phase  raised  the  hydraulic  head  approximate^  0.45  meters. 
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The  injection  occurred  at  a  rate  of  4.85  ni*/d  on  Simulation  Days  15  and  16  at  Row  61,  Cohinm 
1 1,  and  Layer  7.  These  data  were  input  to  the  Well  package. 

F.  SOLUTION  METHOD 

MODFLOW  solves  the  partial  differential  equation  of  groundwatw  flow  using  a  finite- 
difference  technique.  Several  methods  to  solve  the  resulting  system  of  algebraic  equations  are 
available,  including  the  strongly  mq)Hcit  method  (SIP),  sfice  successive  over  relaxation  method 
(SOR),  and  the  preconditioned  conjugate  gradient  method  (PCG).  For  the  present  study,  the 
PCG2  (HiP  1990)  package  was  chosen  because  it  is  efficient  and  because  it  requires 
approximate^  one-fourth  the  storage  needed  by  the  SIP  package. 

The  PCG2  module  solves  the  set  of  linear  difference  equations  iteratively.  The  equations 
are  produced  firom  the  finite-difference  model,  and  can  be  »q)ressed  in  matrix  notation  as 

A  *  X  =  b 

where  A  is  the  coefficient  matrix,  x  is  the  vector  of  hydraufic  heads,  and  b  is  the  vector  of 
defined  flows,  hi  an  iterative  solver,  it  is  assumed  that  the  solution  has  converged  'v^dlen  some 
residual  (difference  in  results  between  successive  iterations)  is  less  than  a  user-^ecified 
convergence  criterion  (Hill,  1990).  For  the  problan  at  hand,  and  for  most  other  problems,  the 
convergence  criterion  was  chosen  to  be  of  the  same  order  of  magnitude  as  the  measurement 
uncertainty.  Hence,  a  value  of  0.01  meters  was  used  since  hydraulic  head  measurements  were 
reported  to  2  deomal  places  (fijr  example,  62.75  meters). 


A  relaxation  parameter  (RELAX)  is  also  q)ecified  as  input  in  the  PCG2  package.  Hill 
(1990)  suggested  that  a  value  of  1.0  to  be  used.  A  smaller  value,  such  as  0.99  or  0.98  may  reduce 
the  number  of  iterations  needed  for  the  solution  to  converge.  Initially,  the  relaxation  was  set  at 
0.98. 
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G.  COMPUTER  SIMULATIONS 


After  all  the  parameters  were  estabKdied  and  the  input  files  were  written  for  each  module, 
MODFLOW  was  run  numerous  times.  The  first  MADE-2  shnulations,  conducted  during  the 
qimmftr  of  1992,  «CTmplified  the  model  to  steady  state  (Gray,  1992).  The  first  transient  simulations 
were  conducted  in  the  summer  of  1993  (Gray,  1993).  Many  of  the  techniques  for  manipulation  of 
data  to  a  usable  form  described  above  were  established  at  that  time.  The  1993  simulations 
included  a  homogeneous  hydraulic  conductivity  field,  based  on  the  findings  of  the  AT-2  pump 
test,  as  well  as  a  more  realistic  heterogeneous  field.  The  early  vertical  discretization  of  the  aquifer 
described  above  refers  to  the  1993  simulations. 

The  present  study  is  a  direct  continuation  of  work  performed  in  the  summer  of  1994  (Gray 
and  Rucker  ,  1994).  Rve  cases  (M2-5-1  through  M2-5-5)  were  run  during  the  sunnner  of  1994 
and  are  described  in  Table  13.  Simulations  M2-5-1  through  M2-5-5  do  not  include  conductivity 
profiles  K-72  throu^  K-81  as  ftiey  were  not  measured  until  ^ring  1995.  Computations  were 
performed  on  a  Sun  Sparc  2  workstatioiL 


TABLE  13.  SUMMARY  OF  MODFLOW  CASES  M2-5-1  THROUGH  M2-5-5. 


Case 

RELAX 

WETDRY 

(metm) 

Specific 

Yield 

Confined 
Storage  Coef. 

Run  Time 
(min) 

Final  Volume 
Error 

1 

0.98 

-0.1 

0.1 

0.0001 

60 

-0.25% 

2 

1.00 

-0.1 

0.1 

0.0001 

NA 

-0.24% 

3 

0.98 

-0.01 

0.1 

0.0001 

72 

-0.25% 

4  ^ 

0.98 

-0.1  1 

0.2 

0.0005 

94 

-1.52% 

5 

0.98 

-0.1 

0.05 

0.00005 

58 

-0.23% 

were  conducted  prior  to  Case  1  (M2-5-1),  however  this  was  the  first  one  to 
converge.  Case  2  tested  the  relaxation  parameter  in  the  PCG2  solver  package.  Thou^  the 
relaxation  parameter  did  little  to  change  the  results,  it  is  befieved  to  have  dowed  convergence. 
Case  3  examined  the  change  in  the  WETDRY  parameter  in  the  BCF2  package.  Tins  parameter 
controls  the  rewetting  of  dry  cells.  A  negative  sign  indicates  the  head  of  the  cell  in  quesrion 
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depends  solefy  on  the  head  of  the  cell  below  it.  The  absohrte  value  of  WETDRY  is  the  excess 
head  required  to  rewet  the  cell,  calculated  as  the  head  in  the  cefl  below  minus  the  bottom 
elevation  of  the  cefl  in  question.  For  example,  if  the  head  of  a  cefl  in  Layer  2  increased  to  a  value 
of  the  absolute  value  of  WETDRY  above  the  bottom  of  an  overlying  dry  cefl  in  Layer  1,  the  cefl 
in  Layer  1  would  rewet.  A  positive  sign  for  WETDRY  requires  that  the  head  in  the  four  adjacent 
cells  in  the  same  layer  must  be  WETDRY  above  their  bottoms  for  rewetting  to  occur.  This 
method  is  unstable  when  there  are  fixed  heads  in  the  grid,  as  demonstrated  by  unsuccessfiil 
emulations. 

Cases  4  and  5  investigated  the  effects  of  uniformly  changing  the  aquifer  storage 
parameters.  Increasing  storage  increased  run  time  and  final  volume  error.  The  volume  error  is  the 
percent  relative  discrepancy  between  inflow  and  outflow  and  is  accumulated  throu^out  the 
simulation.  Errors  below  1%  are  excellent;  acceptable  errors  range  up  to  10%. 

Rgure  16  conq)ares  the  Case  1  head  contours  for  Layer  4  with  the  tq>per  observed  heads 
for  Simulation  Day  270.  Figure  17  gives  the  correqionding  results  for  Layer  9  and  the  lower 
observed  heads.  The  agreement  is  good.  The  pattern  progresses  firom  a  high  gradient  in  the 
southern  portion  of  the  domain  (near  field)  to  a  gradual  dipping  in  the  northeast  (fer  field).  The 
wider  ^amng  in  the  mid  field  is  a  result  of  the  hi^er  conductivity  values. 

To  more  accurately  check  the  results  of  the  simulation,  a  quantitative  approach  was  used 
in  yMck  the  heads  of  a  continuous  monitored  piezometer  were  compared  to  simulated  heads  m 
the  equivalent  location.  Program  WELL<»PH  extracted  the  head  as  a  fimction  of  time  at  a 
sped&c  location  from  the  MODFLOW  output  file.  The  model  boundary  conditions  changed  only 
16  times  in  a  468.day  simulation  so  the  model  cannot  possibfy  predict  the  erratic  day  to  day 
variations  seen  in  the  measured  hydrogr^hs.  Therefore,  for  compaiiscm  the  program 
HYDRCK3tA(ph)  averaged  the  observed  heads  over  each  stress  period  and  conqiared  the  result 
to  the  unaveraged  simulated  heads  (Gray  and  Rucker,  1994).  Figure  18  shows  the  resultant 
hydrograph  of  piezometer  P53A  for  Case  1. 
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Piezometer  P53A 


(ni)  PB3H 
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Figure  18.  Observed  and  Simulation  M2-‘5“l  Hydrographs  of  P53A. 


In  addition  HYDROGRA  calculated  the  minimum,  maximum,  and  root  mean  square 
(RMS)  differences  between  the  time  averaged  observed  head  and  simulated  head.  These 
calculations  are  summarized  in  Table  14  for  Cases  1,  4,  and  5.  Cases  2  and  3  showed  poorer 
results.  Case  5,  with  the  smallest  storage  coefficients,  gave  the  best  agreement  as  indicated  by  the 
RMS  OTor. 

The  remainder  of  this  section  concerns  work  siq)ported  by  the  present  contract.  In  the  M 
of  1995,  Case  5  was  rerun  at  West  Virginia  University  on  a  90  MHz  Pentium  personal  conq)uter 
(PC).  Slight  differences  were  encountered,  eq)ecialfy  in  run  time  and  percent  discrepancy.  The 
PC  took  approximately  one-fourth  the  time  (15  minutes)  to  converge  conq)ared  to  the  Sun 
Sparcstation  2.  However,  the  PC  had  a  -1.40  %  cumulative  volumetric  error  con^ared  to  the 
Sun’s  0.23% .  These  differences  were  considered  inrignificant. 

VTith  the  addition  often  new  conductivity  wells  (K72  through  K81)  in  the  q>ring  of  1995, 
the  transmissivity  and  vertical  leakance  fields  were  re-kriged  with  Geo-EAS  using  the  procedure 
described  previousfy.  These  revised  conductivities  were  the  baas  of  the  M2-8  series  of 
c^ilatintis  Table  15  lists  the  parameters  used  to  create  the  variogram  models  for  the  nine  layers, 
ttirhiHing  model  type,  nugget,  sill,  and  range.  The  table  also  lists  some  m5>ottant  statistical 
parameters,  such  as  mean  and  variance  of  the  natural  logarithms. 


TABLE  14.  OBSERVED  HEADS  MINUS  MODFLOW  HEADS  FOR M2-5-l,4,5. 


Min. 

(m) 

Min. 

(m) 

Min. 

(m) 

Max. 

(m) 

Max. 

(m) 

Max. 

(m) 

RMS 

(m) 

RMS 

(m) 

RMS 

(m) 

Wefi 

Casel 

Case4 

Cases 

Casel 

Cased 

CaseS 

Casel 

Cased 

CaseS 

P53A 

-0.65 

-1.32 

-0.57 

0.74 

0.51 

0.14 

0.329 

0.228 

0.194 

P54A 

-0.53 

-0.84 

-0.37 

0.39 

0.58 

0.30 

0.143 

0.165 

0.136 

PS4B 

-0.42 

-0.78 

-0.17 

0.43 

0.52 

0.43 

0.147 

0.159 

0.143 

P55A 

-0.53 

-0.80 

-0.37 

0.44 

0.50 

0.44 

0.199 

0.204 

0.199 

P55B 

-0.12 

-0.44 

+0.01 

1.01 

1.01 

1.01 

0.374 

0.374 

0.374 

P60A 

-0.51  1 

-0.51 

-1.51 

0.30 

0.38 

0.30 

0.188~^ 

0.188 

0.188 

P61A 

-0.40 

-0.40 

-0.40 

0.36 

0.36 

0.36 

0.188 

0.188 

0.188 

P61B 

-0.39 

-0.39 

-0.39 

0.23 

0.23 

0.23 

0.154 

0.154 

0.154 

-0.44 

-0.69 

-0.35 

0.49 

0.51 

0.40 

0.215 

0.208 

0.197 
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TABLE  15.  KRIGING  PARAMETERS  FOR  TRANSMISSIVITY  AND  LEAKANCE. 


Parameter 

Layer 

Type 

Nugget 

sm 

Range 

Mean 

Trans. 

1 

1.5 

100 

1.43 

5.06 

EISSI^I 

1.0 

4.7 

100 

1.55 

4.93 

1.5 

4.7 

100 

1.84 

5.32 

4 

1.0 

6.5 

200 

1.83 

5.36 

5 

1.0 

5.0 

80 

1.53 

5.13 

6 

nQHjjjlllll^ 

80 

1.24 

4.68 

7 

Exp. 

0.5 

4.0 

120 

1.05 

3.51 

8 

Exp. 

0.5 

2.5 

100 

0.89 

2.35 

9 

0.0 

2.25 

60 

1.76 

2.19 

VCONT 

1 

1.5 

4.0 

100 

1.43 

5.06 

2 

1^^31111 

1.0 

4.3 

100 

1.39 

4.33 

3 

1.5 

4.5 

90 

1.32 

4.89 

4 

1.0 

7.0 

160 

0.96 

5.59 

5 

1.0 

5.0 

70 

0.71 

5.07 

6 

1.0 

4.0 

70 

0.31 

3.75 

7 

Exp. 

0.5 

3.0 

100 

0.18 

2.87 

8 

Exp. 

0.5 

3.0 

90 

0.02 

2.99 

viketQ  Trans.  =  Transmissivity,  Ejq).  =  Exponential,  Sph.  =  Spherical,  and  =  variance 


Figure  19  conq)ares  the  kriged  conductivities  of  layer  3  for  the  cases  of  M2-5  and  M2-8. 
The  difference  in  the  nrid  field  demonstrates  the  in^ortance  of  a  more  exhaustive  data  set  The 
remaining  transmissivity  and  VCONT  contours  for  the  M2-8  simulations  can  be  found  in 
Appendix  B. 
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Simuktion  M2-8-1  used  the  same  parameters  as  M2- 5-1,  except  for  the  aew  set  of 
conductivity  data  in  the  former.  M2-8-1  took  approximately  45  minutes  to  converge  on  a  486  PC 
with  a  speed  of  33  MHz.  The  Layer  4  and  8  head  contours  (Figure  20)  of  day  270  diow  good 
agreement  with  the  observations  in  Figure  14.  The  Layer  4  results  are  quite  similar  to  those  of 
M2- 5-1  in  Figure  16.  However,  the  average  RMS  error  con^ared  to  the  recording  piezometers 
was  0.207  meters  (as  calculated  by  HYDROGRA)  coii5)ared  to  0.215  meters  in  M2-5-1. 
Simulation  M2-8-2  changed  the  storage  coefficients  to  those  used  in  M2-5-5.  The  model  was 
rerun  on  the  486,  converging  after  45  minutes.  The  head  contours  were  virtually  unchanged.  The 
average  RMS  error  for  M2-8-2  was  0.205  meters,  better  than  M2-8-1,  yet  not  as  good  as  M2-5-5 
(0.197  meters).  Snrmlatimi  M2-8-2  was  taken  as  the  new  base  case,  as  it  included  all  updated 
toformation.  Table  16  lists  the  suimnary  of  the  head  differences  for  die  two  cases  of  M2-8-1  and 
M2-8-2. 


TABLE  16.  SUMMARY  STATISTICS  FOR  M2-8. 


ISBH 

Weft 

M2-8-1 

M2-8-2 

M2-8-1 

M2-8-2 

M2-8-1 

M2-8-2 

P53A 

-0.98 

-0.67 

0.35 

0.36 

0.204 

0.198 

P54A 

-0.64 

-0.38 

0.51 

0.51 

0.154 

0.149 

P54B 

-0.47 

-0.47 

0.40 

0.35 

0.177 

0.176 

P55A 

-0.57 

-0.53 

0.45 

0.45 

0.202 

0.202 

P55B 

-0.17 

+0.01 

1.02 

1.02 

0.381 

0.381 

P60A 

-0.51 

-0.51 

0.35 

0.35 

0.192 

0.194 

P61A 

-0.40 

-0.40 

0.36 

0.36 

0.188 

0.188 

P61B 

-0.39 

-0.39 

0.23 

0.23 

0.155 

0.155 

IBBSg 

-0.52 

-0.42 

0.46 

0.45 

0.207 

0.205 
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Other  experimental  amulations,  such  as  those  of  M2-6,  'v^diich  studied  a  two-dimensional 
veraon  of  M2-5,  and  M2-7,  yMah.  studied  the  effects  of  differential  tenqioral  discretization,  were 
no  hei^  in  understanding  the  flow  phenomena  at  the  MADE-2  site.  These  cases  merely  tried  to 
simplify  the  model  for  solving  the  tran^ort  problem,  discussed  later. 

Following  a  suggestion  of  Dr.  C.  Zheng,  simulations  M2-9  tried  an  approach  different  than 
kiiging  to  estimate  the  conductivity  field.  The  polygonal  method,  in  \^dudi  a  polygon  of  mfluence 
is  assigned  to  each  measurement,  produces  a  aep-Eke  conductivity  array.  Program  NGP  (nearea 
grid  point)  teaed  the  conductivity  data  of  the  77  profiles  for  eadi  layer  and  produced  a  blocky 
conductivity  field  in  ^lich  each  node  is  assigned  the  value  of  the  nearea  conductivity 
measurement.  Figure  21  ^ows  a  conductivity  layer  i^duch  is  representative  of  all  the  layers,  since 
each  profile  extends  vertically  downward  in  the  same  horizontal  location.  Simulation  M2-9-1, 
nging  the  same  aorage  coefficients  as  M2- 5-5,  took  approximately  15  minutes  on  a  90  MHz 
Pentium  PC  to  converge.  Figure  22  shows  Layers  4  and  9  on  day  270,  equivalent  to  the  observe 
iqiper  ^md  lower  staged  piezometers  of  March  8,  1991,  in  Figure  14.  Hie  major  disCTepancy 
occurs  in  the  near  field  \^ere  tiie  contours  are  ^read  out  too  mudi,  indicating  a  hi^er 
transmissivity  than  eiqiected.  The  averaged  RMS  error  was  also  hi^er  than  previous  sinaulations 
with  a  value  of 0.232  m.  No  adtHtimial  simulations  were  conducted  using  the  polygonal  method. 

The  output  of  the  flow  model  was  used  as  input  for  the  tran^oit  model,  MT3D.  A 
Knifing  package  called  LKMT18,  written  by  Zheng  (tiie  author  of  MT3D),  produced  an  ou^ut  file 
which  contained  the  heads,  velodties,  and  ceDrby^ceH  flow  terms  needed  to  solve  the  Ihre^ 
.fimpncinnal  tranqiftit  equation.  Hie  next  section  describes  the  process  by  which  the  transport 

equation  was  solved. 
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Figure  2 1 .  Nearest  Grid  Point  Transmissivity  Distribution  for  Layer  9. 


SECTION  VI 

TRANSPORT  MODELING 


A.  THEORY 


The  tran^ort  equation  was  solved  separately  using  MT3D,  a  modular  three-dimensional 
tranq)ort  model  (Zheng,  1990).  Decoupling  the  flow  problem  from  the  tran^ort  problem  is 
permissible  if  density  variations  are  negligible.  However,  if  density  varies  significantly  due  to 
concentration  or  temperature  variations,  the  two  problems  must  be  solved  simultaneous^; 
groundwater  velocities  can  be  affected  by  buoyancy  forces.  The  validity  of  this  assumption  is 
discussed  further  in  a  later  section. 

MT3D  can  ordy  model  only  one  contaminant  at  a  time.  Even  thou^  several  ^edes  were 
injected  in  MADE-2,  they  do  not  react  with  each  oflier.  Therefore,  MT3D  was  used  to  amulate 
the  migration  of  each  contamrnant  ind^endently. 

The  tranqjort  equation  that  MT3D  solves  incorporates  the  terms  representing  advection, 
diversion,  reactions,  and  any  sHiks  and  sources.  Equation  (4)  is  the  partial  differential  equation 
governing  three  itimengional  tran^ort  of  contaminants  in  groundwater 


at 


.  A{v.c)  + 

^  ax, 


(4) 


where;  Xi  refers  to  the  Cartesian  coordinate  axes 
C  is  concentration 
t  is  time 

is  the  hydrodynamic  diversion  (diversion  and  diffiision)  tensor 
Vi  is  the  seepage  velocity 

qs  is  volumetric  flux  of  water  per  unit  volume  of  aquifer  representing  sinks/sources 
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Cs  is  the  contammant  concentration  of  the  water  added  by  anks/sources 

6  is  the  porosity 

is  the  chemical  reaction  term 

1.  Advection 

Advecdon  is  the  migration  of  miscible  contaminants  at  the  groundwater  seq)age  velocity. 
MT3D  solves  the  advective-di^erave-reactive  equation  using  a  mixed  Exilerian-Lagrangian 
scheme.  In  the  Eulerian  approach,  the  tranq)ort  equation  is  solved  using  fixed  nodes,  as  in  the 
finite-diflference  method.  This  method  is  advantageous  in  di^eraon  or  reaction  dominated  flow. 
However,  if  advection  e7q)lains  moa  of  the  contaminant  migration,  the  Eulerian  approadi  is 
susceptible  to  large  numerical  dispersion  and  oscillation  in  the  solution  and  may  requne  small  grid 
q)acing  and  time  aeps  (Zheng,  1990).  Numerical  diq>eraon  causes  a  smearing  of  concentration 
fronts  which  should  have  a  sharp  appearance.  Therefore,  the  Lagrangian  approach  is  often  used 
for  the  advection  dominated  flows  vMsk  exist  in  many  field  conditions.  The  Lagrangian  method 
tracks  moving  particles  and  provides  an  efficient  solution  to  problems  with  sharp  concentration 
fronts. 


The  second  term  in  the  transport  equation  advection  and  can  be 


by  the  T  agrangiflin  approach  in  MT3D  using  eitha  the  method  of  characteristics  (MOC), 
the  modified  method  of  diaracteristics  (MMOC),  or  a  combination  of  the  two  called  the  hybrid 
method  of  diaracteristics  (HMOC).  A  fourth  approach  uses  the  Eulerian  Upstream  Finite 
Differencing  (UD).  Eadi  Lagrangian  method  uses  a  particle  tracking  technique,  \nhich  decreases 
the  amnimt  of  numerical  diversion  and  osdDation  in  flie  sohitioiL 


The  MOC  method  uses  forward  tracking  of  particles,  in  \diich  a  set  number  of  particles 
are  introduced  in  each  computational  cell  of  the  domain  and  are  assigned  a  position  and 
concentration.  The  particles  are  tracked  uang  a  very  small  time  stqi  and  the  average  particle 
concentration  for  the  cell  is  calculated  at  the  end  of  each  time  increment  The  concentrations  of 
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each  particle  are  also  \q)dated  to  reflect  changes  due  to  diversion  and  chemical  reactions.  The 
major  advantage  to  the  MOC  method  is  that  it  is  virtually  free  of  numerical  diversion.  However, 
the  method  is  conqmtationally  difficult  and  requires  large  amounts  of  conqiuter  memory  to  store 
particle  locations.  It  can  also  lead  to  large  mass  balance  discrepancies,  as  it  is  independent  of  the 
princ^le  of  conservation  of  mass. 

The  MMOC  method  also  tracks  particles,  but  does  so  by  placing  only  one  particle  at  the 
node  of  each  cell  anH  tracking  it  backwards  in  time  to  find  its  position  at  the  old  time  level  Since 
only  one  particle  per  cell  is  used  in  the  MMOC  method,  it  is  much  fester  than  the  MOC  method 
whCTe  many  particles  are  sometimes  in  a  cefl.  Moreover,  the  MMOC  method  places  a  new 
particle  at  the  node  in  the  begmning  of  each  new  time  level,  thereby  eliminating  the  need  to  store 
the  location  of  the  particle.  The  method  is  virtually  free  of  numerical  dispersion,  except  at  sharp 
concentration  fronts,  and  is  intended  for  use  ^^ere  sharp  fronts  do  not  occur. 

The  third  method  of  particle  tracking  incorporates  the  advantages  of  both  methods  above. 
Near  sharp  fronts  the  MOC  method  distributes  the  required  number  of  particles  around  each 
front  Elsewhere,  the  MMOC  method  is  employed  to  reduce  the  total  number  of  particles  needed. 
The  solution  automatical^  adapts  to  either  metiiod  as  concentration  fronts  diss^ate  by  di^eraon 
or  chemical  reactions.  A  user-^ecified  criterion  controls  the  switching  between  MOC  and 
MMOC.  In  certain  drcumstances  the  HMOC  method  may  not  give  the  optimal  solution,  in 
case  the  MOC  or  MMOC  may  by  chosen  (Zheng,  1990). 

Lastly,  an  iq>stream  finite  dififerencing  option  for  the  advection  term  is  also  included.  The 
UD  method  may  lead  to  large  numerical  diversion  for  problems  having  sharp  fronts,  but  is  more 
computationally.  The  method  also  conserves  mass,  and  may  reduce  the  mass  balance 

discrq>ancy  at  every  step. 
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2.  Diq)ersion 


Diversion  itt  porous  media  refers  to  the  q>reading  of  contaminants  over  a  greater  region 
than  would  be  predicted  solefy  from  the  Darcy  velocity  vectors.  Hydrodynamic  diversion 
incorporates  both  mechanical  di^asion  due  to  deviations  from  the  average  groundwater  velocity 
within  a  representative  elementary  volume  and  molecular  difriision  due  to  concentration 


d  ^ 

variations.  The  third  term  of  Equation  (4), 

axiV 


represents  diversion  with  being  the 


hydrodynamic  diversion  tensor.  A  fully  e^qifidt,  block-centered  jfinite-difference  method  is  used 
to  approximate  this  term.  The  limitation  associated  with  an  ej^Hcit  solver  is  the  small  time 
stepping  criteria  needed  to  ensure  numerical  stability.  However,  due  to  the  large  memory 
requirement  for  the  particle  tracking,  an  m^licit  solver  would  more  than  likely  exceed  memory 
availability  on  smaller  personal  con^utos. 


3.  GiemicalBeactioiis 

The  fourth  term  in  Equation  (4)  rq)resents  reactions  associated  with  contaminant 
migration  such  as  equilibrium  controlled  linear  or  nonlinear  sorption  and  first  order,  irreversible 
rate  reactions.  The  chemical  reaction  termed  is  e^ressed  as  (Zheng,  1990): 


= 


(5) 


vviiere:  Pb  is  the  bulk  density  of  the  aquifer  material 

C  is  the  concentration  of  the  sorbed  ^edes  per  unit  mass  of  porous  media 
X  is  the  rate  constant  of  the  first  order  rate  reactions 

The  reaction  teim  allows  for  sorption,  the  mass  transfer  process  between  the  contanrinants 
dissolved  in  the  groundwater  and  the  solid  contaminants  sorbed  on  the  porous  media.  It  is 
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assumed  that  the  two-phase  interaction  between  soEd  and  solution  is  in  equflibrium  and  that  the 
reaction  is  6st  enough  to  be  considered  instantaneous.  The  relationshq>  between  the  sorbed  and 
the  dissolved  concaitration  is  caEed  an  isotherm.  MT3D  can  simulate  both  Enear  and  nonEnear 
isotherms.  The  nonlinear  isotherms  available  are  the  Freondfich  and  the  Langmuir. 

Biodegradation  and  radioactive  decay  can  be  simulated  as  first-order  rate  reactions.  Ihe 
rate  constant,  X,  is  usually  given  as  the  half-Efe,  >^ch  is  the  time  needed  to  decrease  the 
concentration  to  one-half  its  initial  value. 

4.  Sinic  and  Sources 

The  fifth  term  of  the  governing  equation,  y  C.,  represents  any  sinks  or  sources  of  sohite 

mass  that  may  leave  or  enter  the  domain.  They  can  be  of  either  point  or  areal  type.  Point  sources 
include  wens,  rivers  and  drains;  areally  distributed  sinks  or  sources  include  recharge  and 

evapotranqnration.  The  concentration  must  be  given  for  any  source  term.  On  the  other  hand,  it 

is  not  necessary  to  specify  the  concentration  for  any  sink,  since  it  is  assumed  that  it  is  equivalent 
to  the  ambient  groundwater  concentration.  The  only  exception  is  evapotran^iration,  in  which 
pure  water  is  taken  from  the  aquifer  (Zheng,  1990).  In  reaEty,  however,  volatile  contaminants, 
mriiiHtng  tritiated  water,  may  leave  through  evapotranspiration.  This  is  a  definite  hmitation  m 
shnulating  MADE-2. 


B.  DISCRETIZATION 

MT3D  is  designed  to  allow  use  of  the  same  grid  as  MODFLOW,  and  this  feature  was 
exploited  in  the  present  study.  Thus  the  MT3D  domain  was  divided  into  66  rows,  21  columns, 
and  9  layers.  The  block  centered  formulation  places  a  node  at  the  center  of  each  cell,  \^ere 
MODFLOW  calculates  head  and  MT3D  calculates  concentration.  The  hydrauEc  and  chermcal 
parameters  such  as  hydraufic  conductivity  and  diq>ersivity  are  assumed  to  be  uniform  over  the 

entire  cell  (Zheng,  1990). 
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The  teii5)oral  discredzatioii  is  also  identical  to  MODFLOW’s  with  the  stress  periods  and 
timft  steps  given  in  Table  10.  However,  the  time  step  used  in  the  in^lidt  solution  of  the  head 
values  by  the  flow  model  may  be  too  large  for  the  e?q>licit  tran^ort  model,  which  places 
restrictions  on  the  time  increments.  Therefore,  the  MODFLOW  time  steps  are  automatically 
divided  by  MT3D  into  transport  steps,  smaller  increments  of  time  in  which  heads  are  kept 
constant  as  the  change  in  concentration  is  calculated. 

C.  INITIAL  AND  BOUNDARY  CONDITIONS 

Because  the  experiment  requires  a  transient  simularion,  initial  conditions  are  necessary  to 
solve  the  transport  equation.  For  the  MADE-2  model,  time  began  14  days  prior  to  any  injection 
of  contaminants.  Hence,  all  concentration  values  at  the  beginning  of  the  tranqiort  shnulation 
were  set  to  zero.  The  boundary  condidons  for  afl  time  were  that  on  the  lateral  boundaries  the 
concentration  of  entering  water  was  set  to  zero  \^ereas  outflowing  water  carried  the 
concentration  of  the  iqistream  node  out  of  the  domain.  It  was  inqificitly  assumed  that  there  is  no 
loss  of  contaminant  throu^  the  water  table,  an  assunqition  vshidi  is  unHkefy  to  be  precisely 
correct. 

D.  TRANSPORT  PARAMETERS 

The  following  aH^triniial  parameters  were  needed  fiir  transport  modehag:  soil  porosity, 
di^ersivity,  diemical  reaction  constants,  and  contaminant  source  strmigths. 

1.  Porosity 

The  porosity  of  the  aquifer  is  needed  to  convert  flie  specific  disdiarge  or  Darcy  velocity 
calculated  by  the  flow  modd  to  seqiage  velodty  for  flie  solution  of  the  tranqiort  equation. 
Porosity  was  measured  for  84  samples  taken  from  only  four  cordioles.  Based  on  these  84 
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an  average  porosity  of  0.32  was  calculated.  Lacking  adequate  data  to  infer  the 
e?q)ected  ^atial  variation  in  porosity,  this  value  was  used  for  every  cell  in  the  domain. 


2.  Di^ersivity 

Boggs  et  aL  (1993)  describe  two  cases  in  which  diq)ersivity  was  calculated  from  the  first 
and  second  q)atial  moments  of  the  observed  MADE-2  tritium  phime,  a  base  case  and  an 
extrapolated  case.  The  base  case  was  calculated  solely  from  the  observations;  the  extrapolated 
case  included  a  model  for  the  portions  of  the  phime  v^iiich  had  q)read  beyond  the  sanqiling 
domain.  Longitudinal  dispersivities  of  19.6  meters  and  9.5  meters  were  calculated  for  the  base 
case  extrapolated  case,  reflectively.  A  horizontal  transverse  difiersivity  of  2.2  meters  was 
calculated  in  only  the  base  case.  This  may  be  an  overestimate  due  to  the  theory  of  the  tranfiort 
model  used  in  calculating  this  number.  Transverse  difiersrvrties  are  usually  taken  to  be 
approximately  10%  of  the  longitudinal  (Freeze  and  derry;  1979).  In  this  woric,  vertical 
transverse  and  lateral  transverse  difiersrvities  were  10%  of  the  longitudinal  value. 

The  molecular  difBision  coeffident  of  tritium  in  water  was  estimated  using  the  Wilke- 

-4  2 

Chang  method.  This  value  was  multiplied  by  an  assumed  tortuosity  of  0.25  to  yield  2. 16x10  m 
/  day,  for  the  molecular  diffusion  coefihdent  of  tritium  in  a  saturated  porous  media.  The  same 
value  was  used  for  the  other  contaminants.  This  carmot  be  justified,  but  its  effect  on  the  results  is 
probably  negligible.  (A  later  recalculation  gave  a  value  of  0.475x10“*  m^  /  day  for  the  molecular 
/tiffiicinn  coefl&doit  of  tririiitn  in  a  saturated  porous  media.  Again,  the  effect  of  the  error  is 

probably  negligible.) 

3.  Chemical  Reactions 

Tritium  in  the  form  of  tritiated  water  is  nonreactive  in  the  groundwater  environment,  but  it 
undergoes  radioactive  decay  with  a  12.26  year  half-life.  Sorption  of  tritiated  water  does  not 
occur,  however  the  hydrocarbon  contaminants  do  ej^erience  sorption.  MT3D  amnlates  the 
effects  of  sorption  in  the  Chemical  Reaction  module  through  the  use  of  a  retardation  fector  (R). 
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A  retardation  fiictor  of  5  inq)]ies  the  contaminant  moves  5  times  slower  than  the  groundwater 
seepage  velocity.  A  linear  equilibrium  isotherm  was  assumed  and  the  retardation  6ctor  was 
calculated  as; 


R 


i.jJ^ 


1  + 


Pb 


Ka 


(6) 


where:  Kd  is  the  distribution  coefficient. 

Pb  is  the  bulk  density  of  the  dry  soil 
Bijjc  is  the  porosity  of  cell  ij  Jc 

Although  Boggs  et  aL  (1993)  obtained  estimates  for  R  from  field  and  laboratory  data,  these  could 
not  be  entered  directfy  to  MT3D.  Instead,  Equation  (6)  was  solved  for  Kj  and  Kj,  0,  and 
Pb  were  entered,  allowing  the  program  to  calculate  R  intemalfy.  Table  17  ^ows  the  values  used 
to  simulate  sorptioiL 


TABLE  17.  PARAMETERS  FOR  ORGANIC  TRANSPORT  SIMULATIONS. 


Organic 

e 

Bulk 

Density 

X  10‘  (*/«3) 

Retardation 

Factor 

Ki  xlO* 

(m^/g) 

RCl 

(day^) 

Initial 

Concentration 

(*/»3) 

Benzene 

0.32 

1.77 

1.30 

5.42 

0.0070 

68.01 

0.32 

1.77 

1.42 

7.59  1 

0.0064 

1  o~BC!B 

0.32^ 

1.77 

1.32 

5.78 

0.0046 

0.32 

1.77 

1.24 

4.33 

0.0107 

41.44  1 

Biodegradation  was  amulated  by  using  experimentally  determined  rate  constants  obtained 
from  MacIntyre  et  aL  (1993).  Degradation  kinetics  were  calculated  from  field  data  and  were 
approximately  first  order.  The  rate  constant  RCl  in  Table  17  was  entered  directly  into  the 

tran^ort  model 
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The  last  cohmm  of  Table  17  lists  the  initial  concentrations  used  for  the  transport 
gifnnIatinTis  These  vahies  were  calculated  by  dividing  the  injected  naass  for  each  organic  by  the 
total  volume  of  solution  (9.7  m^).  The  initial  concentration  is  also  entered  directly  in  the  tran^ort 
model  The  mttial  tritium  concentration  (not  listed)  was  0.0555  Ci/m^ . 

4.  Contaminant  Sources 

The  source  of  contaminants  was  a  line  of  five  injection  wells  ^aced  1  meter  ^art.  In  the 
model  these  were  consolidated  into  one  well  located  at  Row  61,  Column  1 1,  and  Layer  7.  MT3D 
allows  the  user  to  qaedfy  the  concentration  of  injected  contaminants  for  every  stress  period.  For 
the  MADE-2  experiment,  injection  started  on  Simulation  Day  14  and  lasted  until  Day  16, 
corresponding  to  Stress  Period  2.  There  was  no  well  flow  in  any  other  stress  period. 

E.  MT3D  OUTPUT 

MT3D  allows  the  user  to  create  to  five  output  files  including  a  generic  output  file,  an 
unformatted  concentration  file,  a  mass  balance  file,  a  point  observation  file,  and  a  configuration 
file.  The  generic  output  file  contains  all  relevant  information  about  the  simulatioiL  It  reproduces 
the  input  data  as  well  as  the  results.  At  the  end  of  eadi  time  step,  the  generic  output  file 
summarizes  the  mass  budget 

The  unformatted  concoitration  file  and  the  configuration  file  are  created  for  post¬ 
processing  of  the  output  data.  Program  POSIMOD,  created  by  2aieng,  is  used  to  prq)are 
output  in  a  format  for  Surfer  to  contour.  The  unformatted  concentration  file  is  a  binary  file 
containing  only  concentration  values  for  every  layer  at  user  ^edfied  times.  Ihe  configuration  file 
contains  model  discretization  data  needed  by  the  postprocessor  for  graphical  presentation.  The 
generic  oxctprxt,  concentration,  and  configuration  files  are  automatically  created  by  MT3D. 

The  mass  balance  and  observation  files  are  only  created  when  the  option  is  invoked  by  the 
user  in  the  Basic  Transport  module  (BTN).  The  mass  balance  file  lists  the  mass  in,  mass  out. 
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diSerence  between  Tnass  in  atiH  out,  and  relative  percent  diflference  between  the  two  at  every 
tran^ort  stq)  throughoxit  siniulation.  Ute  observation  ffle  lists  the  concentration  at  any  user- 
designated  cell  at  every  tran^ort  step. 

F.  COMPUTER  SIMULATIONS 

The  "ittial  MADE-2  modeling  efforts  of  Gray  (1992)  did  not  include  tranq>ort  modeling. 
Gray  (1993)  used  MT3D  with  the  parameters  given  above  to  simulate  the  tritium  phune  only,  but 
MT3D  ran  so  slowly  compared  to  MODFLOW  that  it  was  not  practical  to  simulate  the  entire 
e:q)eiiment.  FigHty  days  were  simulated  but  only  die  first  52  days  gave  meaningful  resuhs.  For 
the  homogeneous  conductivity  case  the  tranq>ort  emulation  took  approximately  6.4  hours  to  run 
on  the  Sun  Sparc  2.  However,  the  heterogeneous  case  took  almost  three  times  longer  to  simulate 
the  same  time  paiod.  The  homogeneous  case  simnlated  advective  flow  using  the  HMOC  solver, 
vshereas  the  heterogeneous  simulation  used  the  MOC  option.  Gray  concluded  that  heterogeneous 
conductivity  produced  more  realistic  results  than  a  sin^Med  homogeneous  conductivity  solution. 

The  transport  modeling  of  Gray  and  Rucker  (1994),  summarized  in  Table  18,  was  a  little 
more  successful  The  head  solution  from  MODFLOW  case  M2-5.1  was  used  as  input  to  MT3D 
running  on  a  Sun  Sparc  2. 

In  Table  18  the  ‘IKq).’  column  identifies  vtiiich  simnlations  used  di^eraon.  Not  all  runs 
dispersion,  as  it  was  ihou^  to  be  the  cause  of  the  model  not  completing.  The  ‘Long. 
Di^.  (m)’  column  gives  the  longitudinal  di^ersivity  used  in  eadi  simnlation.  Transverse  lateral 
and  transverse  vertical  coii?)onents  of  di^ersivity  were  taken  to  be  one-tenth  this  value.  The 
‘Last  Stm  Day*  the  number  of  the  last  simulation  day  before  the  shnulation  stopped 

because  of  an  overflow  error  or  user  termination.  The  overflow  errors  occurred  in  the  mass 
balance  file.  FORTRAN  has  a  non-zero  miniiinim  and  a  maximmn  number  it  can  interpret,  usually 
10'*^  and  lO'^^S  re^ectively.  When  these  bounds  are  exceeded,  an  underflow  or  overflow  error 
occurs.  For  the  sinmlations  of  Table  18,  extreme  mass  accretion  in  the  mass  balance  file  caused 
most  to  terminate.  Neither  the  mass  balance  calculation  itself  nor  the  cause  of  the 
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accumulation  of  mass  is  fully  understood  by  the  authors.  The  ‘Mass  Discrep.’  column  denotes  the 
discrepancy  between  mass  in  and  out  of  the  domain.  The  ‘Phune  Characteristics  cotumn 
describes  the  general  features  of  the  tritium  phune.  Many  of  the  simnlations  had  some  negative 
concentrations.  TTiis  non  physical  result  is  caused  by  numerical  oscillation  in  the  solution,  and 
occurred  most  frequently  with  the  HMOC  or  UD  methods.  In  general,  the  lower  dispersion 
coefficients  caused  fewer  concentration  values  to  be  negative  and  resulted  in  fester  run  times. 


TABLE  18.  MT3D  SIMULATIONS  BY  GRAY  AND  RUCKER  (1994). 


Run 

Advection 

Method 

Disp. 

Long. 

Disp. 

(m) 

Last 

Sim. 

Day 

Run 

Time 

(hours) 

Mass 

Discrep. 

% 

Plume  Characteristics 

1 

HMOC 

10.0 

30.2* 

15.75 

+7.93 

wide  spread,  some<0 

2 

MMOC 

yes 

10.0 

5.0* 

1.75 

N/A. 

iniection  not  started 

3 

MMOC 

no 

N/A. 

129.4 

10.4 

+82 

wide  spread 

4 

HMOC 

no 

N/A. 

20.4* 

0.72 

+19.2 

not  recorded 

5 

MOC 

no 

N/A. 

62.1* 

3.5 

-13.1 

confined  to  7  cells 

6 

HMOC 

no 

N/A. 

140.9 

17.38 

+17.2 

confined  to  8  cells 

7 

HMOC 

yes 

10.0 

44.6*  1 

47.05 

+4.55  1 

wide  spread,  lots<0 

8 

UD 

yes 

10.0  1 

61.2 

16.6 

—0 

wide  spread,  lots<0 

9 

UD 

yes 

4.0 

90.4 

<21.4  ^ 

~0 

wide  spread,  lots<0 

11 

UD** 

yes 

4.0 

90.4 

<29 

~0 

identical  to  case  9 

12 

UD 

1.0 

128 

<5.37  . 

M) 

realistic.  lots<0 

13 

UD 

yes 

0.0 

138.3 

<8 

realistic,  few<0 

14 

HMOC 

yes 

1.0 

105.9* 

<12.6 

+12.3 

realistic.  lots<0 

♦  run  terminated  by  user  **  double  precision 


The  extension  of  the  summer  of  1994  work  began  wife  Emulations  exploring  different 
dispersion  values  using  a  90  MHz  Pentium  personal  corumter.  However,  no  run  extended 
beyond  Simiilatinn  Day  226.  The  code  was  investigated  feorou^ify  to  find  errors  vdiich  would 
cause  fee  mass  accumulation.  The  exact  source  of  error  was  not  fiiund,  but  it  was  theorized  that 
it  was  related  to  fee  rewetting  capability  of  MODFLOW.  Before  Stress  Period  6  fee  water  table 
was  in  a  steady  decline.  Then  a  large  positive  rediarge  caused  the  water  table  to  rise  abruptly 
approximately  2  meters  in  Stress  Period  7.  In  the  belief  that  the  abrupt  rise  in  fee  water  table  was 
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at  feuh,  a  new  flow  Emulation,  designated  as  M2-6,  divided  the  previous  seventh  stress  period  of 
24  days  into  6  sq>aiate  stress  periods  with  4  days  each  to  give  a  more  gradual  rise  in  head. 
Again,  overflow  errors  caused  the  tranq)ort  simulation  to  end  after  only  219  days. 

fti  order  to  reduce  the  con:q)lexity  of  the  model,  the  grid  was  reduced  to  two  dimensions 
by  nsmg  one  layer  to  represent  the  vertical  discretization.  This  was  the  first  simulation  to  model 
the  entire  468-day  eiqieiiment.  Three  separate  sinmlations  were  run  to  test  the  effects  of 
diqiersivity  and  the  different  advection  solvers.  Table  19  summarizes  the  values  used. 


TABLE  19.  SUMMARY  OF  2-D  TRANSPORT  SIMULATIONS  (M2-6). 


Run 

#  of  Dim. 

Days  Simulated 

Run  Time 

Dl,Dv»Dt  (m) 

Conductivity 
of  Layer 

1 

2 

468 

Ihr 

5,  0.5,  0.5 

1 

2 

2  ' 

468 

1  hr 

10,1,1 

1 

3 

2 

468  1 

1  hr 

10, 1,  1 

4 

Figures  23  and  24  diow  Runs  1  and  2  at  Simulation  Days  148  and  344,  re^ectively.  All 
concentration  values  were  normalized  by  the  injected  concentration.  Tliese  runs  used  the 
hydraufic  conductivity  values  of  Layer  1  of  the  nine-layer,  three-dimensional  grid.  Tliese  results 
do  not  compare  well  with  the  observed  contours  of  Rgure  25.  Die  longitudinal  spread  of  the 
phune  is  not  as  fir  as  the  observed,  only  going  to  a  maximum  of  75  meters  on  Day  468 
(not  shown).  In  addition,  the  lateral  spread  is  excessive.  The  conductivity  was  thou^t  to  be  the 
hindrance.  Therefore  Run  3  used  flie  conductivity  of  Layer  4.  Figure  26  shows  the  Run  3 
normalized  tririiitii  concentration  on  Days  148  and  344.  Diere  is  very  little  m5)rovemenL 
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27  days  132  days  224.  days  328  days 


Figure  25.  Observed  Tritium  Relative  Concentration 
at  59.5  meters  MSL.  Somce:  Boggs  et  al  (1993). 


Day  344. 


Dr.  Thftfig  the  author  of  MT3D,  suggested  reducing  the  value  of  the  Percel  variable,  read 
from  the  Advection  module  input  file.  Percel  is  the  Courant  number,  the  maximum  number  of 
cells  any  particle  is  allowed  to  move  in  one  tran^ort  step.  The  Courant  number  estabh^es  the 
ttiavimnni  tfme  stq)  allowed  in  Order  to  maintain  stability  in  the  numerical  procedure.  In  the 
manual  for  MT3D,  it  is  suggested  that  Percel  lie  between  0.5  and  1.0.  However,  Dr.  Zheng 
suggested  ngmg  about  0.01.  The  first  successfiil  three-dimensional  tranq)ort  simulation,  using  a 
Percel  of  0.01,  took  approximatefy  7  days  on  a  90  MHz  Pentmm  PC.  It  quit  on  Day  446  (17 
stress  periods),  due  to  an  error  in  an  input  file  for  the  last  stress  period.  The  simulation  used  the 
UD  method  for  solving  the  advection  term  and  had  a  longitudinal  diq)eraon  coefficient  of  5.0 
meters.  The  flow  gntmlatinn  of  M2-5-5  was  used  as  input  for  the  seepage  velocities.  Figure  27 
diows  the  normalized  tritium  contour  of  Simulation  Day  344.  This  coireq)onds  to  the  observed 
day  of  328  since  the  sanulation  started  14  days  prior  to  the  injection.  Ahhou^  the  simulated 
plume  did  not  migrate  as  frr  as  tihe  observed  plume  in  Figure  25,  it  looked  much  more  reasonable 
than  previous  attenq)ts.  More  simulations  were  performed  to  increase  the  migration,  as 
summarized  in  Table  20. 


TABLE  20.  SUCCESSFUL  HORIZONTALLY  ISOTROPIC  TRANSPORT  SIMULATIONS. 


Run 

#  of 

Dim. 

Days 

Simulated 

Run 

Time 

Advection 

Solver 

Dl  y  Dt  y  Dv 

(m) 

Po'cd 

DTO 

0 

1 

3 

446 

7  days 

UD 

5,0.5  ,0.5 

0.01 

0.001 

0.32 

2 

3 

226 

13  days 

MOC 

10.0, 0. 1,0.1 

0.01 

0.001 

0.32 

3 

3 

302 

7  days 

MOC 

2.0,  0.2, 0.2 

0.1 

0.01 

0.29 

H 

3 

468 

4  days 

UD 

2.0, 0.2, 0.2 

0.1 

0.01 

0.29 

86 
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Runs  2  and  3  used  the  MOC  advection  solver,  but  did  not  produce  desirable  results. 
Many  concentration  values  were  less  than  zero,  indicating  large  numerical  oscillations.  The  MOC 
solver  produced  high  numerical  diversion  in  the  fiir  field  at  sharp  concentration  fronts.  Plots 
after  Snmilatinn  Day  148  show  increastng  instability  as  the  phirne  migrates  ferther  downstream. 
The  MOC  option  also  took  much  longer  than  the  UD  to  simulate  a  given  period.  Only  runs  using 
the  UD  option  generated  concentration  values  that  were  always  positive.  The  DTO  value, 
q)ecified  by  the  user,  is  the  largest  tranq>ort  step  the  simulation  can  use.  As  the  di^ersivity 
decreased,  so  did  the  Courant  number  and  the  maximum  transport  step  to  ensure  stability. 

The  last  column  of  Table  20  lists  the  porosity  (0)  used  in  the  simulation.  As  porosity 
decreases,  seepage  velocity  increases.  The  porosity,  while  remaining  within  the  range  of  the 
TTiAaqirp-niflnts  was  decreased  in  the  last  two  simulations  to  try  to  pudi  the  plume  beyond  the  75 
meter  distance  obtained  with  the  original  porosity.  Though  not  shown  here,  the  normalized 
tritium  phime  of  Run  4,  with  a  porosity  of  0.29,  did  not  migrate  as  fax  as  the  previous  simulation 
due  to  the  decreased  dispersion.  However,  it  did  produce  a  skinnier  plume,  much  more  realistic 
than  its  predecessor. 

All  above  «ni.ilatinn.s  used  the  flow  model  of  M2-5-5.  TTie  M2-8  flow  simiilations,  wiiich 
included  the  latest  conductivity  data,  were  conducted  after  the  above  tranqiort  emulations. 


The  gmmiatiftn  of  the  trrtmm  phune  was  also  attenqited  using  the  M2-9-1  flow  model, 
based  on  the  nearest  grid  point  hydraulic  conductivity  fidd.  The  shnulatitm  lasted  approximate^ 
4  days  on  a  90  MHz  Pentium  and  used  a  longitudinal  disperavity  of  1.0  meters.  Tlie  upstream 
finite-differencing  method  was  used  to  solve  the  advection  term,  and  the  Courant  number  was  set 
to  0.01.  Figure  28  diows  the  normalized  tritium  phune  for  Day  344.  The  phune  onfy  readied 
approximately  25  meters  downgradient  from  the  source.  Obviously,  the  contaminant  was  stuck  m 
a  very  low  conductivity  block  and  did  not  readi  areas  of  hi^er  conductivity. 
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SEcnoNvn 

EFFECTS  OF  HORIZONTAL  ANISOTROPY 


As  described  in  the  previous  section,  several  problems  hampered  attempts  to  simulate  the 
MADE-2  tritium  phime.  First,  MT3D  Version  1.80  ran  excessively  slowly,  up  to  3  weeks  in 
some  qmnlarinns  This  was  personalty  firustrating  as  wdl  as  conqiutationalty  expensive.  Second, 
the  MT3D  particle  tracking  methods  produced  exces^e  numaical  diversion.  Evai  with  low  to 
zero  hydrodynamic  diq)ersi<m,  the  concentration  of  cells  located  6r  downstream  often  increased 
above  the  injected  concentration  with  the  MOC  method.  Last,  the  simulations  Med  to  accurately 
match  the  observed  plume  for  runs  that  did  conqilete  the  468-day  experiment.  The  simulated 
phune  did  not  q)read  6r  enough  in  the  longitudinal  direction.  The  dtyping  of  the  calculated  heads 
to  the  northwest  in  the  fer  field  also  seemed  to  contradict  the  observed  plume  behavior,  which 
migrated  to  the  northeastern  comer  of  the  domain.  All  these  issues  had  to  be  dealt  with  if  a 
simulation  was  to  be  completed. 

A-  MT3D  CODING  MODMCATIONS 

The  excessive  run  time  of  the  tran^ort  code  MT3D  seemed  related  to  the  rewetting 
capability  of  the  MODFLOW  flow  code.  Before  the  Courant  number  was  pushed  to  a  very  low 
munber,  decreasing  the  time  stepping  fector,  MT3D  usually  stopped  around  Day  226.  However, 
iniTMiigrir.  mass  accTction  Started  around  Day  162  vHam  an  influx  of  recharge  started  to  raise  the 
water  table  and  rewet  dry  cells  in  the  tqiper  layers.  The  subdiviaon  of  original  Stress  Period  7 
into  shorter  stress  periods  in  order  to  smooth  the  abrupt  increase  in  head,  Med  to  enable  the  code 

to  run  any  longer. 

Koch  (1994)  reported  finding  a  dimension  error,  wiiidi  essentially  acted  as  a  virus,  in  the 
advection  package  of  MT3D.  But  after  moditying  MT3D,  Version  1.1,  his  shnulations  stffl  only 
ran  to  about  140  days,  crashing  as  a  result  of  overflow  errors.  Further  investigation  by  Kodi 
(1994)  identified  an  additional  problem;  the  artificialty  small  time  step  criterion  for  the  sink/source 
terms  at  the  water  table  (Equations  4.29  throu^  4.32  of  Zheng,  1990).  The  MODFLOW  flow 
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model  was  set  up  in  such  a  way  that  prior  to  Day  130,  a  negative  recharge  (indicating 
evaporation)  was  used.  After  this  time,  an  abnq)t  addition  of  positive  recharge  lead  to  the  rise  of 
the  water  table,  rewetting  previous^  dry  cells.  This  is  wiiai  the  problem  of  small  time  steps 
occurred.  Koch  (1994)  described  the  difficulty,  and  fixed  the  problem  for  the  particle  tracking 
MOC,  MMOC,  and  HMOC  options.  Koch’s  modified  MT3D,  referred  to  as  Version  1.1,  was 
ftsg^tial  in  the  following  simulations.  Even  though  Koch’s  corrections  allowed  MT3D  to  run  in 
a  reasonable  rime  frame,  his  rririum  plumes  only  extended  to  approximately  75  meters  after  224 
Simulation  Days,  con^ared  to  the  observed  225  metms  (Koch,  1994). 

Zheng  (1995)  produced  results  in  which  the  tritium  phime  extended  to  about  120  meters 
downgradient.  However,  his  amulations  were  for  an  “averaged”  steady-state  condition  in  >^ch 
the  rime  stepping-problem  did  not  occur.  Modeling  only  the  tritium  plume,  he  tested  various 
diqiCTsivity  values  as  well  as  different  conductivity  variations.  He  concluded  that  the  simulated 
plume  is  more  senative  to  the  way  the  hydraulic  conductivity  field  is  genoated  than  to  the 
di^ershdty  values  used,  and  recommended  using  the  nearest  grid  point  method  instead  of  kriging 
the  conductivities  in  order  to  reduce  smoothing. 

B.  HORIZONTAL  ANISOTROPY 

All  previous  modeling  attempts  had  Med  to  accuratefy  match  tiie  observed  phime.  The 
plume  never  extended  fiuther  tiian  120  meters  downgradient,  i^ereas  the  observed 
pinmft  extoided  to  about  225  meters.  Recalling  the  results  of  pun^  tests  AT-1  and  AT-2, 
horizontal  anisotropy  was  ^ipfied  to  the  conductivity  field  to  increase  longitudinal  advecbon. 


1.  Row  Modeling  with  Horizontal  Anisotropy 

MODFLOW  allows  the  user  to  specify  horizontal  anisotropy  in  the  Block-Centered  Flow 
package  as  long  as  the  princqial  axes  of  hydraulic  conductivity  are  aligned  with  the  coordinate 
axes  and  the  anisotropy  is  constant  in  each  layer.  The  horizontal  anisotropy  ftctor,  TRPY,  is  the 
ratio  of  column  transmissivity  (or  hydraufic  conductivity)  to  row  transmissivity.  One  value  of 
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TRPY  may  be  given  for  each  layer.  A  vahie  of  1.0  produces  horizontally  isotropic  conditions.  A 

value  greater  <’ban  1.0  increases  the  column  transmissivities  vritile  leaving  the  row  transmissivities 
unchanged. 

The  horizontally  anisotropic  simulations  were  based  on  the  conductivity  field  of  the  M2-8 
simulations  \^di  incorporated  the  new  conductivity  readings  fi:om  wells  K72  -  K81.  Several 
values  of  TRPY  were  tested  in  order  to  increase  longitudinal  velocities  and  plume  advection. 
Table  21  lists  the  results  of  the  uniformly  horizontally  anisotropic  flow  simnlations. 


TABLE  21.  UNIFORMLY  HORIZONTALLY  ANISOTROPIC  FLOW  SIMULATIONS. 


Simulation 

M2-8-2 

M2-8-3 

M2-»4 

M2-8-5 

M2-8-6 

TRPY 

1.0 

2.0 

10.0 

50.0 

100.0 

RMSofPSSAtml 

0.1985 

0.1991 

0.2036 

0.2082 

0.2112 

0.1487 

0.1427 

0.1377 

0.1457 

0.1518 

RMSofP54B(m) 

0.1760 

0.1752 

0.1841 

0.2096 

0.2016 

RMSofPSSAtm) 

0.2030 

0.1993 

0.1998 

0.2091 

0.2160 

0.3815 

0.3685 

0.3514 

0.3567 

0.3671 

0.1943 

0.1906 

0.1849 

0.1846 

0.1869 

0.1883 

0.1883 

0.1883  1 

0.1883 

0.1883 

0.1550 

0.1550 

0.1550 

0.1550 

0.1550 

Average  RAfS  (m) 

0.2056 

0.2023 

0.2000 

0.2071 

0.2097 

%  Discrepancy 

3.84 

3.70 

2.90 

1.40 

0.84 

Five  uniformly  amsotrqpic  flow  amulations  were  conducted,  with  anisotropy  (TRPY) 
ranging  from  1  to  100.  lie  last  two  values  of  50  and  100  were  used  only  to  observe  the  head- 
gradioit  changes  noticed  from  tiie  hi^  Darcy  velocities.  To  date,  no  study  has  beai  found  to 
tndi^?*p!  that  such  extreme  conditions  exist  anywhere  in  the  world;  therefore  no  physical  meaning 
can  be  attached  to  the  hi^er  numbers  (Rucker,  1995).  Ahhou^  the  highest  horizontal 
anisotropic  value  found  in  the  literature  was  only  3.5,  at  a  site  in  Dawsonvffle,  Georgia  (Madia 
and  Randolph,  1987),  a  model  value  of  10.0  gave  the  best  averaged  RMS  error.  The  RMS  error 
is  the  root  mean  square  difference  between  the  averaged  observed  head  and  the  simulated  head. 
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These  flow  shnulatioiis  were  conducted  on  a  33  MHz  486  PC,  and  took  approxinaatefy  45  minutes 
to  con^ilete.  Figures  29  through  33  show  heads  for  simulations  M2-8-2  throu^  M2-8-6  for 
Layers  4  and  9  at  Day  270,  corre^onding  to  the  iq>per  and  lower  kriged  observed  heads  of 
March  8,  1991,  which  are  shown  in  Figure  14.  The  head  contours  look  reasonabfy  accurate 
except  for  the  later  gimiilatimis  of  M2-8-5  and  M2-8-6  which  had  unreaBsticalty  high  anisotropy. 
The  most  noticeable  problem  is  with  the  contour  of  64.20  meters  MSL.  The  conductivities  for 
this  region  show  high  variability  through  a  lateral  cross-section. 
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Throu^  careful  examinatioii  of  individual  RMS  error  values  for  each  simulation,  it  was 
noticed  that  every  piezometer  had  a  dififerent  “optimum”  TRPY  value.  For  exan^le,  piezometer 
P53A  had  a  gmaller  difference  in  head  vsiien  anisotropy  was  equal  to  1,  vvhereas  piezometer  P60A 
was  optimized  at  a  value  of  10.  This  prompted  an  addidmial  simulation  with  anisotropy  varying, 
not  only  by  layer,  but  also  by  column  and  row.  hCnor  adjustments  were  made  to  the  on^nal 
MODFLOW  source  code  to  read  the  TRPY  value  from  a  two-dimensional  array  reader,  with  one 
array  for  each  layer.  This  allows  each  ceE  to  have  a  different  TRPY  value.  The  modified  version 
of  MODFLOW  was  saved  as  MF95.  MF95  was  demonstrated  to  be  correct  by  reproducing 
previous  gHrmlatinTis  The  coding  dianges  can  be  found  in  v^endix  C. 

The  new  sinnilation,  M2-8-7,  used  a  smooth^  varying  anisotropy  fector  (TRPY)  array  for 
each  of  the  nine  layers,  produced  by  kriging  optimum  values  of  anisotropy  at  the  eight  q)ecified 
piezometer  locations,  The  kriged  values  ranged  from  1.0  to  50.0,  and  the  resulting  TRPY 
isopleths  can  be  seen  in  Figure  34.  An  i^per  kriged  TRPY,  taken  from  piezometers  designated  by 
“A”  was  assigned  to  Layers  1  through  4.  The  lower  kriged  TRPY,  taken  from  piezometers 
designated  as  “B”  was  assigned  to  Layers  5  throng  9.  The  flow  simulation  took  approximately 
45  to  compile  on  the  486  PC.  Figure  35  shows  the  heads  of  Layers  4  and  9  at 

Simiilatinn  Day  270.  The  contours  compare  well  with  the  observed  heads  in  Figure  14. 

M2-8-7  produced  the  best  averaged  RMS  error  to  date,  with  a  sHghtfy  lower  value  of 
0.1944  meters.  However,  the  real  advance  was  the  increased  velocity  fidd.  Program  VELCOMP 
compared  the  Darcy  velocity  arrays  of  the  isotrtpic  emulation  (M2-8-2)  to  that  of  flie  variably 
anisotropic  gimiilatian  (M2-8-7)  and  computed  a  percent  relative  increase  at  all  nodes.  An 
average  increase  of  23 1%  was  calculated  fbr  flie  longitudinal  velodties. 
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Figure  34. 


2.  Tra3iq)ort  Modeling  of  M2-8-7 


MT3D  Version  1.1,  modified  by  Koch,  was  used  to  solve  the  tranq>ort  equation  for  the 
variably  horizontally  anisotropic  flow  field  of  M2- 8-7.  An  additional  change  suggested  by  Koch 
to  q)eed  the  simulations  was  to  sk^  the  tran^ort  calculations  during  the  first  10  days  of  the  flow 
dttniiatiftn  because  there  were  no  contaminants  present  before  Day  14.  The  first  simulations 
considered  the  tritium  phune.  Virtually  afl  parameters  remained  the  same  as  before,  including 
poroaty,  diffiisivity,  and  decay  coefl&dent.  The  onfy  major  diange  was  the  increase  of  the  mitial 
tritium  concentration  from  0.0555  Ci/m^  to  1000  Ci/m  to  avoid  underflow  errors.  As  long  as  all 
terms  used  in  calculating  the  change  in  concentrations  are  linear,  the  above  increase  is  acceptable, 
because  the  calculated  concentrations  are  divided  by  the  mitial  concentration  to  obtain  a 
normalized  tritium  concentration.  A  longitudinal  diq)eravity  of  0.5  meters  md  transverse  and 
vertical  dispersivities  of  0.01  meters  were  used  on  the  advice  of  Koch.  The  MOC  option  was 
used  to  solve  the  advective  portion  of  the  transport  eq[uation  because  Koch  (1994)  reported  that 
the  Ollier  particle  trackers  gave  rise  to  numerical  diversion  and  smaller  time  stepping  at  sharp 
concentration  fironts,  and  the  upstream  finite-difference  option  had  not  been  corrected. 


The  first  tran^ort  simulation  in  this  series  took  ^proximatefy  7  hours  to  conq>lete  on  a 
Sun  Sparcstation  10.  Figure  36  diows  the  contours  for  normalized  tritium  concentration  at 
approximately  100-day  intervals.  Hie  exact  days  were  diosen  to  match  the  snapshots.  The  last 
day.  Day  456,  does  not  have  a  corresponding  sn^shot,  but  is  shown  to  demonstrate  the  migratian 
of  tritium  at  the  end  of  the  ej^eriment.  It  was  obvious  that  the  predicted  tritium  phime  was  still  a 
poor  fit  to  the  observed  data  shown  in  Figure  25  because  the  center  of  mass  migrated  away  from 
the  injection  zone,  fa  reality  the  contaminaiit  remained  primarily  in  the  low  conductivity  zone 
near  the  injection  site  vMe  a  tail  of  low  concentrations  extended  to  the  nrid  and  6r  field  regions. 
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C.  ANISOTROPY  VARYING  BY  CONDUCnvnY 


Based  on  geologic  considerations.  Young  (1994)  had  proposed  a  geological  history  of  the 
area  around  the  MADE-2  site,  which  posited  braided  streams  giving  rise  to  a  metastable 
meandering  river  channel  Boggs  et  aL  (1992)  had  concluded  that  the  MADE-2  site  contained 
lenses  of  extreme  heterogaieity.  The  deposit  characteristics  of  a  meandering  river  as  described  by 
Freeze  and  Cherry  (1979)  almost  identically  match  the  characteristics  of  the  MADE-2  ^e, 
including  coarse  sand  and  gravel  along  point  bars  and  large  amounts  of  clay  and  sflt  from  channel 
fill  deposits.  From  these  considerations  Young  (1994)  concluded  that  a  charmel  may  have  existed 
through  the  site  and  that  sediments  were  deposited  in  such  a  way  that  p«meability  is  hi^er  in  the 
intigTfiiHtnal  direction  of  the  channel  than  in  the  lateral  directiotL  Figure  5  show^  an  aerial 
photograph  of  the  site  with  the  channel  axis  cutting  through  the  site  at  about  30  east  of  north. 

With  inciglit  gained  from  the  previous  set  of  calculations,  a  new  MODFLOW  simulation, 
M2-8-8,  was  performed.  Rehfeldt  et  aL  (1992)  had  suggested  that  the  presence  of  the 
hypotheazed  channel  in  the  mid  field  may  be  in  the  regions  of  IQi  >  10  cm/s  (~8.6  m/d). 
Tharefore,  to  emulate  the  effects  of  the  channel,  a  high  anisotropy  value  was  assigned  to  the 
higher  conductivity  cells,  and  low  anisotropy  fiictor  to  lower  conductivity  cells  for  run  M2-8-8. 
Program  TRPY  tested  the  conductivity  arrays  fi)r  each  layer  and  wrote  the  new  anisotropy  arrays 
to  the  Block  Centered  Row  package  input  file.  After  several  tests,  it  was  concluded  that  cells 
v^ose  horizontal  conductivity  was  above  3.0  m/d  (—3.5x10  csa/s),  ^ould  be  have  TRPY  set  to 
15,  and  cells  with  lower  conductivity  should  have  TRPY  equal  to  1.  Figure  37  shows  two  typical 
layers  with  the  high  anisotropy  cells  shaded  in  gray. 

The  flow  model  was  rerun  to  obtain  new  head  and  Darcy  velocity  values  fi^r  the  tran^ort 
gtmniatifin  The  Calculated  averaged  RMS  head  discrqiancy  for  the  simulation  was  0.1958 
meters,  slightly  higher  than  M2-8-7,  yet  still  better  than  simulation  M2-5-5  of  Gray  and  Rucker 
(1994).  The  predicted  head  contours  of  Figure  38  are  a  reasonable  fiicsimile  to  the  observed  ones 
in  Figure  14. 
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Again,  trithim  was  the  first  contannnant  to  be  tested  for  solute  ndgratioiL  All  other 
parameters  described  above  were  adopted  for  this  siinulation  as  well  The  resuhant  phime 
contours  of  Figure  39  show  very  good  agreement  with  the  observed  (Figure  25).  The 
contaminant  moves  slowly  in  the  near  field,  increasing  ^eed  as  it  alters  the  mid  regions  of  the 
domain.  The  last  days  show  a  long  tail  trailing  off  to  the  northeast,  just  as  observed. 
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Figure  37.  Regions  of  TRPY  =  15  Indicated  By  Shading.  Left:  Layer  3;  Right:  Layer 
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The  txan^ort  of  the  four  organics  were  also  tested  with  the  new  flow  model  All 
parameters  of  the  above  simulation  were  used  for  these  gmulations  with  the  exception  of  the 
sorption  and  biodegradation  coefficients.  Table  17  diows  the  values  used  to  calculate  the 
retardation  fectors  and  the  initial  concentrations  for  each  of  the  four  organic  contaminants.  The 
table  does  not  include  any  information  about  This  contaminant  was  not  modeled,  since  much 
of  the  information  needed  for  simulations  was  not  found  in  any  of  the  literature.  The  coefficient 
of  molecular  Hiffiisirm  in  a  saturated  porous  medium  was  set  to  the  value  for  tritium  for  aU  of  the 
organics.  This  is  incorrect  but  the  resulting  errors  are  probably  small 

Each  dissolved  organic  tran^ort  simulation  took  between  4  to  7  hours  to  complete  on  a 
Sun  Sparcstation  10.  Figures  40  through  44  show  the  normalized  concentrations  on  a  vertical 
profile  through  the  longitudinal  midsection  of  the  domain  for  Days  44,  148,  240,  344,  and  456.  A 
YZ  profile  through  column  11  was  used  to  present  the  contours  to  conD?)are  to  the  snapshot  data 
of  Boggs  et  al  (1993)  presented  in  Figures  45  through  49. 

Again,  there  is  good  agreement  between  observed  and  simulated  plumes.  The  migration  of 
each  amulated  phime  seems  to  follow  the  path  of  its  corre^onding  observed  plume.  Early  in  the 
gfnrmlatinn  on  Day  44,  the  contaminants  stay  in  the  near  field;  however  the  sinnilated  plumes  are 
more  concentrated  in  the  lower  portion  of  the  domain.  Since  the  injection  took  place  in  Layer  7 
(57.5  meters  MSL)  and  the  water  table  decreased  for  the  first  130  days  of  simulation,  it  seems 
logical  that  the  amulated  contaminant  does  not  rise  to  around  60  meters  MSL,  as  compared  to 
the  observed  phime.  The  poative  net  recharge  in  the  seventh  stress  period  (Day  148),  i^hich 
raised  the  water  table,  also  puHed  the  bulk  of  the  contaminant  to  the  upper  portion  of  the  domain. 
This  affected  the  tritium  plume  the  most,  allowing  a  better  fit  to  the  observed  tnrium  phime  of 
Figure  46.  The  positive  rediarge  fiiroughout  mudi  of  the  remainder  of  the  simulation  kept  the 
majority  of  eadi  phune  aroimd  60  meters  MSL. 

The  ctmniatftH  longitudinal  ^read  of  each  contaminant  did  not  match  the  observations 
perfectly,  but  did  reproduce  the  gross  patterns  of  tranqiort.  The  simulated  tnrium  phime,  for 
example,  migrated  to  about  105  met^  downgradient  firomihe  injection  point  on  Day  148  (Figure 
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41).  Its  observed  cotinterpart  only  migrated  to  about  65  meters  (Figure  46).  The  next  aiapshot 
shows  the  observed  tritium  phune  (Figure  47)  to  q)read  out  to  about  200  meters,  with  the 
simulated  phime  moving  only  out  to  150  meters.  This  discrepancy  probably  results  ftom  the 
hydraulic  conductivity  field  used  in  simulating  the  flow  and  tran^ort  of  the  chemicals.  Although 
kriging  allows  the  modeler  to  acquire  the  best  estimate  of  the  variability,  without  having  a 
conqiletely  exhaustive  data  set  an  exact  match  is  iiDpossible. 

The  dispersivity  values  in  the  transport  calculations  are  the  only  other  major  source  of 
uncertainty  in  the  modeling.  Dispersivity  accounts  for  spreading  due  to  heterogeneity  ^^*ich 
causes  variations  in  the  flow  velocities  and  paths.  However,  the  di^ersivities  assumed  for  these 
simulations  seemed  adequate.  Just  as  in  Zheng  (1994),  it  is  concluded  that  the  simulated  plumes 
are  more  sensitive  to  the  treatment  of  hydraulic  conductivity  than  to  the  di^eravity  values  used. 

Figures  50  through  54  present  areal  displays  of  the  dissolved  organic  plumes  for  Layer  4 
for  Simulation  Days  44,  148,  240,  344,  and  456.  These  figures  do  not  have  corresponding 
observation  plots,  but  they  diow  how  the  contaminants  are  tran^orted. 
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Nonnalized  Tntium  Plume  of  Column  1 1 ,  Day  44  Normalized  Benzene  Plume  of  Column  1 1 ,  Day  44 
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Figure  40.  M2-8-8  Normalized  Concentration  Profiles  Along  Column  1 1 ,  Day  44. 
Isopleths  Presented  at  0. 1,  0.01 , 0.001,  and  0.0001. 


Normalized  Tritium  Plume  of  Column  1 1 ,  Day  1 48 
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Figure  41.  M2-8-8  Normalized  Concentration  Profiles  Along  Column  1 1,  Day  148 
Isopleths  Presented  at  0. 1,  0.01,  0.00 1,  and  0.0001. 


Normalized  Tntmm  Plume  of  Column  1 1 ,  Day  240  Normalized  Benzene  Plume  of  Column  1 1 ,  Day  240 
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Figure  42.  M2-8-8  Normalized  Concentration  Profiles  Along  Column  1 1 ,  Day  240. 
Isopletbs  Presented  at  0.1, 0.01, 0.001,  and  0.0001. 


Normalized  Tritium  Plume  of  Column  1 1,  Day  344 
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Normalized  Tritium  Plume  of  Column  11 ,  Day  456 
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Figure  44.  M2-8-8  Normalized  Concentration  Profiles  Along  Column  1 1 ,  Day  456. 
Isopleths  Presented  at  0.1, 0.01, 0.001,  and  0.0001. 
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C/Co 

Observed  Normalized  Concentration  Profiles,  Day  44.  Source;  Boggs  et  al.  (1993) 
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Figure  46,  Observed  Normalized  Concentration  Profiles,  Day  148.  Source;  Boggs  et  al.  (1993). 


Tritium  Benzene 


^  CN  O  CD  CD  ^ 

ID  CO  (D  tn  in  lO 


^  CM  O  CD  CO  ^ 
CO  CD  CO  kH  CO  lO 


(uj)  uoipAai3  (uj)  uojpAaj] 


118 


Figure  47.  Observed  Normalized  Conceutration  Profiles,  Day  240.  Source:  Boggs  et  al.  (1993). 
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Figure  48.  Observed  Normalized  Concentration  Profiles,  Day  344.  Source:  Boggs  et  al.  (1993). 
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Figure  49.  Observed  Normalized  Cooceutration  Profiles,  Day  456.  Source:  Boggs  et  al.  (1993). 


BENZENE  NAPHTHALENE  p-XYLENE 


121 


Figure  50.  M2-8-8  Normalized  Organic  Concentrations  in  Layer  4,  Day  44. 
Isopleths  Presented  at  0.0001. 
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Figure  5 1 .  M2-8-8  Normalized  Organic  Concentrations  in  Layer  4,  Day  148. 
Isopleths  Presented  at  0.001  and  0.0001. 
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Figure  52.  M2-8.8  Normalized  Organic  Concentrations  in  Layer  4,  Day  240 
Isopleths  Presented  at  0.01, 0.001,  and  0.0001. 
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Figure  53.  M2-8-8  Normalized  Organic  Concentrations  in  Layer  4,  Day  344. 
Isopleths  Presented  at  0. 1, 0.01, 0.001,  and  0.0001. 
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Figure  54.  M2-8-8  Normalized  Organic  Concentrations  in  Layer  4,  Day  456. 
Isopleths  Presented  at  0.01, 0.001,  and  0.0001. 


Program  MASSCALC  was  written  to  calculate  the  zeroth,  first,  and  second  spatial 
moments  of  the  concentration,  and  to  find  the  total  mass,  center  of  mass,  and  variance  of  the 
contaminant  plumes.  Corre^onding  moment  calculations  on  fidd  data  (Boggs  et  aL,  1993)  used 
three-dimensional  numerical  integration  to  account  for  the  variation  of  concentration  between 
san^ling  points.  However,  in  a  finite-difference  calculation,  the  concentration  is  assumed 
constant  over  the  entire  cell  Therefore,  MASSCALC  calculated  the  moments  by  summation. 
Tables  22(a)  through  22(e)  list  the  resuks  of  MASSCALC  fi:om  simulation  M2-8-8. 


TABLE  22  (A).  TRITTIIM  PLUME  CHARACTERISTICS  FROM  SIMULATION  M2-8-8. 


Day: 

44 

148 

240 

344 

456 

MTMo 

.9931 

.9982 

.9818 

1.443 

.9306 

Xc(m) 

1.06 

1.97 

4.88 

9.367 

13.25 

Yc(m) 

1.61 

6.17 

27.06 

93.38 

159.51 

Zc(m) 

57.12 

58.36 

59.42 

59.61 

58.85 

Gxx^(m^) 

4.73 

7.61 

23.90 

68.94 

178.25 

Ovy^  (m^) 

16.93 

92.61 

1584 

4081 

3775 

oj  (m") 

.579 

1.09 

1.33 

1.57 

1.62 

TABLE  22  (B).  BENZENE  PLUME  CHARACTERISTICS  FROM  SIMULATION  M2-8-8. 


Day: 

44 

148 

240 

344 

456 

M/Mo 

1.02 

ATI 

.257 

.150 

.0369 

Xc(m) 

1.11 

1.64 

2.22 

7.04 

6.62 

Yc(m) 

0.9 

3.95 

8.26 

30.6 

64.1 

Zc(m) 

57.09 

57.81 

59.60 

60.77 

59.44 

Oxjc^  (m^) 

4.69 

6.34 

9.19 

15.1 

21.4 

(m^) 

8.15 

30.6 

219 

1490 

4170 

yy 

cj  (m") 

.45 

.852 

1.15 

1.24 

.499 

126 


TABLE  22  (C).  NAPHTHALENE  CHARACTERISTICS  FROM  SIMULATION  M2-8-8. 


44 

148 

240 

344 

456 

M/Mo 

.977 

.502 

.309 

.158 

.0422 

Xc(m) 

1.09 

1.50 

2.79 

4.56 

3.67 

Ye(m) 

.867 

3.19 

4.93 

13.9 

7.64 

Zc(m) 

57.04 

57.84 

58.84 

59.43 

59.53 

(m^) 

4.57 

5.61 

6.91 

7.49 

4.85 

Ow  (m") 

5.88 

9.12 

12.02  1 

331.04~^ 

10.43 

a J  (m^) 

.411 

.651 

1.19 

1.27 

.335 

TABLE  22  (D).  P-XYLENE  PLUME  CHARACTERISTICS  FROM  SIMULATION  M2-8-8. 


Day; 

44 

148 

240 

344 

456 

M/Mo 

.78 

.264 

.095 

.036 

.006 

Xc(m) 

1.17 

1.75 

3.37 

5.18 

4.11 

Ye(m) 

1.03 

3.82 

8.00 

20.76 

6.09 

Zc(m) 

57.11 

57.97 

59.10 

59.76 

59.64 

(m^) 

4.86 

6.46 

7.43 

10.70 

3.65 

a^"(m") 

9.10 

25.56 

133.3 

564.8  ^ 

4.27 

cj  (m"  ) 

.4767 

.859 

1.31 

1.29 

.229 

TABLE  22(E).  O-DCB  H.UME  CHARACTERISTICS  FROM  SIMULATION  M2-8-8. 


Day: 

44 

148 

240 

344 

456 

M/Mo 

1.12 

.648 

.1378 

.0528 

.0087 

Xc(m) 

1.17 

1.63 

3.37 

5.18 

4.11 

Yc(m) 

1.03 

3.85 

BUB 

20.76 

6.09 

Zc(m) 

57.11 

57.91 

59.10 

59.76 

59.64 

(m") 

4.86 

6.25 

7.43 

10.70  ' 

3.65 

Ow^(m") 

9.10 

26.54 

133.3 

564.8 

4.27 

oj  (m") 

.4767 

.799 

1.31 

1.29 

.229 

M/Mo  is  total  mass  to  injected  mass,  Xc  Ye,  and  Zc  are  the  plumes’  centroids  in  the  re^ecrive 
directions,  and  Oxx*,  Oyy*.  and  are  the  piuM^al  second  moments. 
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Tbe  results  in  Table  22  generalfy  conq>are  reasonably  well  to  the  observed  plume 
parameters  from  Boggs  et  al  (1993),  seen  in  Table  23. 


TABLE  23  (A).  OBSERVED  TRITIUM  PLUME  CHARACTERISTICS. 


TABLE  23  (B).  OBSERVED  BENZENE  PLUME  CHARACTERISTICS. 


Day: 

44 

148 

240 

344 

456 

M/Mo 

0.92 

0.43 

0.23 

0.07 

0.06 

Xc(m) 

-0.2 

-1.1 

-0.8 

-1.0 

-1.1 

Yc(m) 

3.8 

6.3 

12.4 

7.7 

7.9 

Zc(ni) 

58.13 

58.68 

59.24 

58.90 

58.67 

Ow"  (m^) 

9.2 

38.4 

826 

24.7 

21.4 

axx^(m^) 

7.9 

6.9 

6.1 

8.4 

10.7 

oj  (m") 

1.9 

1.1 

1.2 

1.2 

1.3 

TABLE  23  (C).  OBSERVED  NAPHTHALENE  PLUME  CHARACTERISTICS. 
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TABLE  23  (D).  OBSERVED  P-XYLENE  PLUME  CHARACTERISTICS. 


Day. 

44 

148 

240 

344 

456 

M/Mo 

0.77 

0.58 

0.23 

0.03 

0.01 

Xc(m) 

-.2 

-1.0 

-1.2 

-1.1 

Yc(m) 

3.6 

6.1 

6.4 

10.5 

WBM 

Zc(m) 

58.04 

58.60 

58.93 

58.66 

57.93 

9.1 

22.8 

16.3 

213 

11.8 

oj  (m") 

7.2 

5.8 

6.4 

4.6 

4.0 

1.9 

1.1 

1.3 

1.0 

1.3 

TABLE  23  (E).  OBSERVED  O-DCB  PLUME  CHARACTERISTICS. 


Day 

44 

148 

240 

344 

456 

M/Mo 

0.81 

0.80 

0.60 

0.21 

0.13 

Xc(m) 

0.0 

-1.6 

0.0 

-0.2 

-1.5 

Yc(m) 

3.9 

7.4 

34.7 

22.1 

8.1 

Zc(m) 

58.09 

58.57 

58.93 

58.90 

58.68 

Cvv^  (m*) 

11.3 

64.4 

3540 

1180 

21.9 

Oxx*  (m*) 

7.9 

7.8 

25.2 

22.4 

7.1 

2.1 

1.4 

2.0 

1.6 

1.8 

The  idative  «»«««  bdancc  (M/Ii^)  for  escii  contamnunt,  coiuj>oted  MASSCALC,  was 
plotted  with,  the  observed  idative  mass  balance  versos  Shnnlation  Day  in  Figures  55  thzoush  59. 
AH  diowed  good  agreement  between  the  two  with  the  exception  of  tiitiam.  Since  Uitium  does 
not  undergo  substantial  decay,  the  mass  balance  should  be  approximately  1.0  (or  100%). 
However,  Figme  55  diows  errors  in  both  observation  and  smuladon.  The  observations  gave  a 
152%  overestimate  at  Day  44  whidi  Boggs  et  aL  (1993)  attributed  to  preferential  sampling  in 
relativefy  hi^  conductivity  zones.  The  23%  underestimate  at  Day  344  of  the  observed  data  was 
due  to  truncation  of  the  leading  edge  of  tire  plume.  The  simulated  tritium  plume  gave  excdlent 
results  except  for  a  44%  overestimate  on  Day  344.  This  simnlation  was  repeated  many  times; 
however,  the  number  persisted  hr  every  calculation.  It  is  not  known  what  caused  tins  error.  A 
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potential  cause  is  the  use  of  the  MOC  option,  which  does  not  inherentfy  conserve  mass.  But  in 
view  of  the  anomalous  nature  of  this  discrepancy,  this  explanation  seems  unlikely. 
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Mass  Balance  (M/Mo) 


Relative  Mass  Balance  for  lyitium 


1  51  101  151  201  251  301  351  401  451 

Snce  Injection 


Figure  55.  Observed  and  M2-8-8  Mass  Ratio  for  Tritium. 
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Relative  Mass  Balance  for  Benzene 


51  101  151  201  251  301  351  401  451 

Days  Since  Injection 


I^gure  56.  Observed  and  M2-8-8  Mass  Ratio  for  Benzene. 
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Figure  58.  Observed  and  M2>8-8  Mass  Ratio  for  p-X\ie3ie. 


Mass  Balance  (M/Mo) 
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Figure  59.  Observed  and  M2-8-8  Mass  Ratio  for  o-DCB. 


135 


FmaHy,  the  center  of  mass  trajectories  for  each  contaminant  were  plotted  to  ^ow 
differences  between  the  observed  and  simulated  plumes.  Figures  60  throu^  64  diow  these 
trajectories  projected  on  the  XY  plane.  In  reality,  the  center  of  mass  moves  vertically  between  57 
and  6 1  meters  MSL.  The  numbers  on  the  plots  represent  the  following  code: 

1  =  injection  period,  simulation  day  16 

2  =  Emulation  day  44 

3  =  simulation  day  148 

4  =  Emulation  day  240 

5  =  Emulation  day  344 

6  —  simiilatinn  day  456 

Eadh  simulated  center  of  mass  is  general^  close  to  the  observed,  ahhou^  the  Emulated 
plumes  deviate  systematically  to  the  east.  The  effects  of  biodegradation  are  seen  in  the  iqpgradient 
motion  of  the  center  of  mass  in  the  later  simulated  plumes  of  n^hthalene,  p-x^ime,  and  o-DCB. 
Similar  behavior  was  noted  in  the  later  observed  plumes  of  benzene,  p>x^ene,  and  o-DCB.  The 
reasons  for  the  lack  of  upgradient  center  of  mass  motion  in  the  predicted  benzene  and  observed 
naphthalene  plumes  are  unclear.  The  case  of  benzene  is  particulaify  puzahng  because  of  the 
excellent  mass  balance  agreement  shown  in  Eigure  56. 

The  emulated  center  of  mass  trajectories  veer  east,  as  opposed  to  the  observed 
trajectories,  t\hich  deviate  sH^tl^  to  the  west  The  maxinmm  discrq>ancy  can  be  seen  in  tritium 
with  an  qyproximate  10  meter  shifi  to  the  east. 
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Center  of  Mass  Trajectory  for  Naphthalene 


-50  0  50 

X(m) 

Figure  62.  Observed  and  M2-8-8  Center  of  Mass  Trajectory  for  Naphthalene. 


D.  FURTEIER  CONSIDERATION  OF  ANISOTROPY 


Although  the  use  of  horizontal  anisotropy  varying  with  conductivity  gave  phunes  \\4ich 
con^)ared  well  with  observations,  the  justification  for  this  procedure  remained  weak.  The 
assumed  Tnavinmni  principal  axis  was  in  the  Y  direction,  12°  west  of  north.  This  does  not  agree 
with  the  resuhs  of  AT-2  (35  °  west  of  north)  or  with  die  hypothesized  channel  axis  (30  °  east  of 
north).  Indeed,  there  is  every  reason  to  befieve  that  the  piinc^al  directions  will  be  dififerent  in  the 
channel  zone  than  elsev\here. 

Another  question  related  to  the  magnitude  of  the  anisotropy.  The  horizontal  anisotropy 
foctor  used  to  increase  velocities  in  the  longitudinal  direction  was  5.8  times  higher  than  estimated 
firom  pump  test  AT-2.  This  discrepancy  was  a  major  concern,  and  additional  research  was 
conducted  in  order  to  justify  the  value  used. 

TTamtiifih  (1966)  describes  the  effects  of  a  fiilfy  penetrating  well  in  an  infinite  nonlealQ^ 
horizontally  anisotropic  aquifer.  The  equal  drawdown  contours  are  ellipses,  with  their  major  axes 
oriented  along  the  maYinnini  prin^>al  direction  of  conductivity.  Furthermore,  he  concludes  that 
the  methods  for  obtaining  aquifer  parameters  for  isotropic  conditions  are  not  valid  for 
anisotropic  conditions.  The  effective  isotropic  values  are  rdated  to  the  actual  prhu^al  values  as 
follows: 


Transmissivity  =  T  = 

(7) 

4iit  1 - 

Storativify  =  S  = 

(8) 

where  T^  and  Ty  are  the  princpal  transmissivities  hi  the  X  and  Y  directions,  re^ectivdy 

u  is  equal  to  - 

4Tt 

r  is  the  radius  to  the  observation  well 
t  is  drawdown  time 
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The  effective  isotropic  transmissivity  is  the  square  root  of  the  product  of  the  two  princ^al 
transmisavides.  Appfying  these  piinc^les  to  the  flow  model  is  a  difficult  task,  since  the  bordiole 
flowmeter  hydraulic  conductivity  values  are  somehow  averaged  in  all  horizontal  directions. 
MODFLOW  uses  the  conductivity  values  to  calculate  the  conductance  betweoi  adjacent  cells.  In 
MODFLOW  the  conductance  is  the  product  of  hydraulic  conductivity  and  cross  sectional  area  of 
flow  divided  by  the  distance  between  the  adjacent  nodes.  The  conductance  equation  incorporates 
anisotropy  by  mult^lying  the  column  transmissivity  by  TRPY  to  get  the  row  transmissivity. 
When  a  TRFY  value  greater  than  one  was  entered  in  the  previous  amulations,  the  overall 
effective  conductivity  was  increased. 


In  an  attenq>t  to  correct  this  problem,  MODFLOW  was  modified  to  maintain  the  observed 
effective  conductivity  by  automatically  reducing  the  conductivity  in  flie  X  direction  to  conqiensate 
for  the  increased  conductivity  along  the  Y  directioiL  The  column  conductance  was  mult^lied  by 
the  square  root  of  TRPY  and  the  row  conductance  was  divided  by  the  square  root  of  TRPY. 
Assuming  CR{^i^^  is  the  conductance  between  nodes  iJ4c  and  ij+14c  and  CC^i/zj^  is  the 
conductance  betweoi  nodes  y,k  and  iri-l  j,k  leads  to: 


CR  ,  2  DELCi  X 


CC  .  2  DELRj  X 


•nku’TRy.i 


u 


TR.  ^  DELR^i  +  TR;^ 


ijtu  J^ELRj 


TC,i,  DELC^,  +  TC^y^  DELCj 


(9) 

(10) 


\^drere  DOJRj  and  DELCi  ate  the  row  and  column  widths  of  cell  ij,k,  respectively 

TRg^  and  TCij^  are  the  row  and  column  transmLsstvities  of  cell  y,k ,  re^ectively 

If  the  domain  is  divided  into  a  square  cells,  as  in  the  MADE-2  modd,  then  DELR  =  DELC,  and 
the  above  two  equations  reduce  to  the  harmonic  mean  of  the  transmissivities  between  two 
adjacent  cells: 
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(11) 


Siibstitutmg  Equation  (13)  into  (16)  fields: 

(17) 

(18) 
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1 


TC;;v  =  T:;t  X 


(19) 


Similarly,  the  invmrse  of  Equation  (13)  can  be  substituted  into  Equation  (16)  yielding: 


^  TRPY 

(20) 

II 

(21) 

Tiix  = 

^TRPY 

(22) 

=  T;,^  Vtrpy 

(23) 

The  transmissivities  of  Equations  (23)  and  (19)  are  substituted  into  Equations  (11)  and  (12)  to 
obtain  the  new  conductances. 


CR  , 


2  X 


ij+^4; 
2 


VTRPY  X  VTRPY 

VTRPY  +  VlRPY 


CR  1  = 


ij+-^ 
2 


cc  ,  = 


2  X  VTRPY  X 

2  X  ^  VTRPY  ^  VTRPY 

T.  1  +T  1 

^  Vtrpy  Vtrpy 


cc  1  =  2  X 


VT^ 


(24) 

(25) 

(26) 

(27) 


Equations  (25)  and  (27)  were  incorporated  into  MODFLOW  and  several  values  of  TRPY 
were  ^iplied  to  the  flow  model  to  test  the  new  changes.  The  first  run  tested  the  old  value  of  15 
widi  the  .same  conditions  as  M2-8-8,  ^^hidi  worked  quite  w^  in  simulating  the  plumes.  The  flow 
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smiulation  took  ^roximatefy  15  mmutes  to  coiiq)lete  on  a  90  MHz  Pentium  PC.  The  heads 
matched  &iily  well  to  the  observed,  with  an  average  RMS  of  0.2127  meters.  However,  the 
simulation  of  tridum  migration  was  unsatis&ctory.  On  Simulation  Day  344  the  tritium  plume 
(Figure  65)  onfy  reached  75  met^s,  mudi  less  dian  the  observed  225  meters  in  Figure  25.  Hie 
reason  for  the  ^ortcoming  was  due  to  the  reduced  Y  conductivities.  A  foctor  of  15  was  used  in 
the  original  M2-8-8  amulations,  but  this  emulation  used  the  square  root  of  15  to  calculate  the 
conductances  between  rows  (CR).  This  obvious^  did  not  produce  the  longitudinal  vdoddes 
needed  to  pu^  the  phime  outward.  Also,  the  effect  of  reducing  the  conductance  between 
columns  (CC)  by  the  square  root  of  TRPY  had  virtually  no  infiuence  on  the  longitudinal 
^reading.  It  onfy  kq»t  the  plume  from  ^reading  in  the  lateral  direcdtm. 

Other  TRPY  values  were  applied  to  the  flow  model  to  increase  veloddes  to  better  match 
the  observed  tritium  plume.  Table  24  lists  the  corrected  transmisavity  simuladons.  The  square 
root  of  the  TRPY  is  provided  to  ^ow  the  actual  value  that  was  muh^hed  and  divdded  into 
observed  values  to  give  the  row  and  column  conductances,  respective^. 


TABLE  24.  SIMULATIONS  OF  M2-8-9  (CORRECTED  ANISOTROPIC 

TRANSMISSIVITY). 


Run 

TRPY 

Vtrpy 

Ava*^ed 

RMS(m) 

Longitudinal  Spread  of  Tritium  Hume  on 
Simulation  Day  344 

1 

15 

3.87 

0.2127 

75  m 

2 

75 

8.66 

0.2216 

100  m 

3 

150 

12.24 

0.2235 

130  m 

4 

225 

15 

0.2243 

crashed 

5 

300 

17.32 

0.2245 

crashed 
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The  flow  models  for  the  M2-8-9  simulations  became  progressively  worse  as  TRPY 
increased.  The  higher  TRPY  Actors  inoreased  the  longitudinal  velocities,  yet  decreased  the  lateral 
velocities  enou^  to  skew  the  flow  directions.  Figures  66  through  70  show  the  heads  of  each  run 
for  Layer  4  and  Day  298,  centered  on  survey  date  April  4, 1991  (Figure  71).  The  velocity  vectors 
are  overlain  on  the  plot  to  ^ow  the  flow  direction.  The  vectors  were  calculated  by  the  program 
VECTOR,  whidi  extracted  the  seepage  velocity  vectors  in  the  X  and  Y  directions  fi’om  the 
MODFLOW  output  flle  and  confuted  a  magnitude  and  direction  for  every  cefl.  The  velooties 
are  actually  coitq>uted  in  MODFLOW  at  each  cell  6ce,  but  were  plotted  at  the  cell  center.  Figure 
72  shows  the  same  layer  and  day  for  simulation  M2-8-8  for  comparison  with  ITgures  66  through 
70.  Even  though  the  amulations  of  M2-8-9  preserve  the  same  hydraulic  gradient  over  the  entire 
domain  as  M2-8-8,  the  individual  contours  do  not  match.  For  example,  the  isopleth  representing 
the  hydraulic  head  of  63.60  meters  for  each  M2-8-9  emulation  extends  into  the  &r  field,  around 
200  meters  downgradient  firom  the  injection.  The  coiTeq)onding  contour  for  M2-8-8  remains  in 
the  near  field.  The  ^pe  of  the  contour  was  also  a  major  concern.  The  hi^er  TRPY  foctors 
produced  ‘humps”  in  the  isopleths.  These  are  obvious  artifocts  due  to  the  treatmoot  of 
transmisavity  and  are  not  seen  in  die  kriged  observed  heads  of  Figure  71  (^lil  4, 1991). 
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Heads  of  Layer  4,  Day  298,  Run  4  of  M2-8-9 


Figure  69.  Head  and  Seq)age  Velocity  in  Layer  4,  Day  298,  M2-8-9,  Run  4. 
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Heads  of  Layer  4,  Day  298,  Run  M2-8-8 
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The  transport  model  was  also  run  for  all  5  shnuktions  of  M2-8-9  to  test  the  tritium 
migration  against  the  successM  M2-8-8  dmulations.  The  last  cohmm  of  Table  24  indicates  how 
fax  the  phune  extended  on  Day  344.  This  value  is  based  on  the  longitudinal  spread  of  the 
normalized  contom  of  0.0001.  In  general  as  the  TRPY  kctor  increased,  the  phune  migrated 
kither  downgradient.  However,  there  was  a  limit  to  the  value  of  the  anisotropy  kctor  vddch 
could  be  used.  The  fourth  run  was  the  first  to  see  the  ill  effects  of  an  increaang  anisotr<q>y.  The 
shniilafinn  became  very  imstable  in  the  earl^  days  and  by  Day  344  (Figure  73)  the  errors 
propagated  throughout  the  entire  domain.  Litense  investigation  of  the  output  file  created  by 
MT3D  ^ows  the  ^ors  starting  around  Day  148,  \\Mch  is  in  the  seventh  stress  period.  Note 
from  Table  12  that  the  transirinn  from  negative  to  positive  net  recharge  began  in  Stress  Poiod  7. 
The  combination  of  skewed  conductivities  with  the  reversal  of  net  recharge  may  have  rendaed 
the  solution  unstable. 

The  fifth  run  of  sinmlation  M2-8-9  also  crashed  around  Day  148.  However,  errors  of  the 
last  two  runs  began  in  two  sq>arate  locations  of  the  domain.  Run  4  started  to  corrtq)t  in  Column 
21,  Row  58,  A^ereas  Run  5  started  in  Column  1,  Row  57.  Both  cells  began  increasing  in 
concentration  frster  than  thdr  nei^bors,  essentially  creating  mass  from  nothing.  Both  cells  were 
also  on  the  boundary.  Hantiish  (1966)  warns  against  applying  the  inverse  of  Equation  (7)  to  some 
flow  fidds.  The  transformation  to  an  isotropic  media  may  produce  an  erroneous  solution  vdien 
the  transformation  of  conductivity  adverse^  affects  the  statement  of  boundary  conditions 
(Hantush,  1966).  Further  research  is  needed  in  order  to  prqpe^  model  tiie  effects  of  horizontal 
anisotropy. 
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SECTION  vm 
DISCUSSIONS 


A.  BUOYANCY  EFFECTS 

Istok  and  Huu^lirey  (1995)  investigated  the  possibility  of  buoyancy-induced  flow  at  small 
concentrations  in  a  two-well  lab  experiment.  Their  study  was  the  first  to  observe  fiiat  relative^ 
low  doisity  differences  (Ap/p  =  7.5  x  10'^)  may  cause  significant  phune  sinking.  The  test  used 
bromide  (Br")  with  a  wide  range  of  injected  concentrations  (50,  250,  500,  750,  and  1000  mg/L). 
Die  centroid  of  each  phune  sank  due  to  negative  buoyancy.  The  implications  may  have  great 
significance  fi>r  numerical  modeling,  especially  for  the  MADE-2  project  ^ce  the  software  used 
assumes  the  denies  of  the  tracer  solution  and  the  ambient  groundwater  are  equal 

The  e?q>eiiment  by  Istok  and  Ifomphrey  (1995)  was  conducted  in  a  steady-state, 
homogeneous  flow  field.  The  test  conditions  were  designed  to  ensure  horizontal  flow  with  no 
vertical  gradients.  The  effects  of  buoyancy-induced  flow  can  be  seat  more  easily  in  this 
e7q)erimait  as  ideal  conditions  were  maintained.  For  the  MADE-2  experiment  the  effects  are  less 
obvious.  The  large-scale  heterogeneities  combined  with  file  fluctuating  water  table  may  nia.sk  the 
buoyancy  effects. 

The  relative  densities  for  the  five  contaminants  were  calculated  (Table  25),  but  none 
separate)^  fell  within  file  range  of  investigation  by  Istok  and  Bfemphrey  (1995).  Their  combined 
total  of  1.5  X  10*^  does  exceed  the  7.5  x  10'^  relative  density  found  to  cause  sinking  Calculations 
using  a  coiqiled  flow  and  transport  model  are  needed  to  determine  \^h.ether  buoyancy  effects  were 
significant  in  MADE-2.  MODFLOW  and  MT3D  are  incapable  of  doing  sudh  calculations,  but 
programs  like  HST3D  and  SUTRA  are. 


TABLE  25.  RELATIVE  DENSITIES  OF  CONTAMINANTS. 
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Contaminant 

Relative  Density  (Ap/p) 

Benzene 

6.8  X  10'^ 

Naphthalene 

7.2  X  10^ 

p-X)iaie 

4.1x10'^ 

0-DCB 

3.3  X  10'^ 

Tritium 

3.8  X  10'^^ 

Total 

1.5  X  lO"^ 

B.  CaOD  CONVERC^CE 

A  grid  refinement  stutfy  was  attenq>ted  in  order  to  determine  if  the  MODFLOW  results 
were  indq>endent  of  horizontal  cell  aze.  Based  on  the  general  theory  of  Ridiardson 
extrapolation,  grid  refinement  studies  are  generally  accepted  in  the  CFD  (Conq)iitational  Fluid 
Dynamics)  community  as  a  tool  to  verify  the  accuracy  and  solution  order  of  a  numerical  solution. 
However,  Roa^e  (1993)  states  that  the  theory  bdiind  the  method  can  be  independoit  of  the 
equations  being  used  and  the  dimenaonalify  of  the  problem.  One  can  easily  ^ly  the  concept  to 
any  numerical  model  as  a  postprocessor  to  solutions  on  two  grids  with  no  reference  to  die  codes, 
algorithms,  or  governing  equations  produced  the  solutions  (Roache,  1993).  Therefi>re  it 
A/^^s  proposed  to  stiufy  the  groundwater  flow  solution  obtained  by  the  MODFLOW  code  with  two 
or  more  grid  sizes. 

Roadie  (1993)  proposes  a  grid  convergence  index  (GCI)  to  account  fi»  the  solution’s 
gridding  as  well  as  its  order.  The  idea  behind  tiie  GCI  is  to  qq)roximatefy  relate  the  error 
obtained  by  an  arbitrary  grid  refinement  study  to  the  error  that  would  be  obtained  by  a  grid 
refinement  study  of  the  same  problem  with  the  same  base  grid  by  a  second  order  method  (p  =  2) 
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and  grid  doubling  (r  =  2).  The  grid  refinement  ratio,  r,  is  the  ratio  of  the  coarse  grid  spacing  (h2 ) 
to  the  fine  grid  pacing  (hi):  r  =  hs/hi . 


The  GCI  was  derived  by  calculating  the  estimated  finctional  error  for  the  fine  grid  solution 
(fi),  that  would  q)proximately  have  the  same  error  with  p  =  2  and  r  =  2.  The  fractional  error  m 
Equation  (29)  is  e7q)ressed  as  Ei. 


(29) 


w&ere  s  = 


(30) 


and  £2  is  the  coarse  grid  solution.  The  GCI  for  the  fine  grid  solution  is  expressed  as 

Sisl 


GCI,  = 


(r'-l) 


(31) 


Notice  that  i^en  the  grid  is  doubled  Ibr  a  second  order  problem,  the  GCI  =  |  e  |. 


The  error  estimator  of  Equation  (29)  does  not  indicate  whether  a  good  estimate  of  the 
solution  has  bear  achieved.  The  error  usua%  sought  is  one  tiiat  gives  confidence  in  the  solution, 
or  a  bound  on  the  error  to  showifthe  solution  is  too  high  or  too  low.  Ei  does  not  provide  sudi  a 
confidence  interval  However,  based  on  cumulative  expmence  in  the  CFD  community,  at  least  a 
marginal  confidoice  level  exists  for  the  error  of  Equation  (30),  obtained  by  a  second  order 
accurate  solution  and  grid  doubling.  Therefore,  GCI  relates  any  grid  refinement  stu%  to  one  of  p 
=  2  andr  =  2. 


Equation  (31)  may  be  evaluated  using  point-by-point  values  or  solution  fimctionals,  such 
as  the  integrated  Darcy  velocity.  Howev^,  the  grid  sizes  must  lie  witiiin  the  asymptotic  range  (as 
assumed  in  the  Ta^or  series  expansion  on  i^ch  the  Richardson  extr^olation  is  based).  This  can 
easily  be  tested  in  a  foiify  straightforward  manner,  as  long  as  the  solution  order,  p,  is  uniform. 

If  the  exact  solution  model  problem  is  known,  tiiea  p  can  be  monitored  by 
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Ep  =  error  /  tf 


(32) 


as  h  is  refined.  The  consistency  of  £p  wiU  verify  the  order  and  indicate  that  the  asymptotic  range 
has  been  adiieved.  Since  this  is  rarely  the  case,  three  grid  solutions  must  be  performed  and  two 
GCl  calculated  firom  the  solutions.  Equation  33  equates  the  GCI  of  the  coarse-intermediate 
calculation  to  the  fine-intermediate  calculation  and  indicates  that  the  asymptotic  range  has  been 
achieved. 

GCIj3  «  r'’GCIi2  (33) 

\^ere  GCI23  is  file  grid  convergence  index  between  an  intermediate  grid  solution  and  a  coarse 
solution 

GCI12  is  the  grid  conv^ence  index  between  a  fine  grid  solution  and  an  intermediate 
solution 

The  technique  was  apphsd  to  several  MODFLOW  emulations  by  using  diSerent  grid 
pacings.  The  standard  flow  shnulations  used  a  grid  pacing  of  5  x  5  meter  cdls  in  the  XY  plane. 
The  solution  of  a  finer  grid  u^g  a  3  x  3  meter  ceE  would  not  converge,  so  onfy  coarser  grid 
pacings  were  used.  Additional  shnulations  were  conducted  by  increasing  the  grid  size  to  a  10  x 
10  and  15  x  IS  meter  ceDs.  Vertical  discretization  was  not  coarsened;  the  MODFLOW  layers 
were  based  on  stratigrphy  of  the  site  and  were  lefi  intact  for  the  grid  refinement  stu^  (as 
suggested  in  conversation  with  Dr.  Roadie).  To  test  the  GCI  of  each  simulation,  file  heads  were 
extracted  for  piezometers  P54A,  P55A,  and  P55B  at  day  270.  Table  26  lists  file  heads  firom  each 
grii 


TABLE  26.  EXTRACTED  HEADS  FOR  VARIOUS  OOD  SIZES. 


WeU 

P54A 

P55A 

P55B 

Ax=  5  m 

64.285  m 

64.204  m 

64.379  m 

0.20561 

Ax  =  10  m 

64.34  m 

64.36  m 

64.36  m 

0.20817 

Ax=  15  m 

64.35  m 

^.38  m 

64.39  m 

0.23136 
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Hie  averaged  RMS  differeaces  of  heads  were  also  tested  for  overall  grid  convergence. 
The  GCls  were  calculated  in  EXCEL  and  Table  27  lists  the  results. 


TABLE  27.  GCI  VALUES. 


1 

P54A 

P55A 

PSSB 

0.00342 

0.0097 

0.00118 

0.0498 

0.00037 

0.0007 

0.00112 

0.2673 

Only  piezometer  P55B  lies  in  the  a^oi^totic  range.  The  other  GCI  calculations  show  a 
difference  of  at  least  1000%  (ten  times  greater),  thus  indicating  that  the  solutions  do  not  lie  in  the 
a^mq)totic  range.  Therefore  the  GCI  calculations  are  meaningless  at  these  sites.  The  GCI 
calculations  for  PS  SB  on  the  other  hand  indicates  a  0.118%  dif^oice  between  the  fine  and 
intermediate  grids  and  a  0.112%  difference  betweoi  the  intermediate  and  coarse  grids.  This 
iiq>lies  that  the  difference  in  head  due  to  using  the  different  grid  pacings  is  negligible  at  die 
location  of  piezometer  PSSB.  The  GCI  was  only  calculated  for  head  differences,  therefore 
nothing  can  be  inferred  about  the  convergence  of  other  calculated  parametos  sudh  as  velocities. 

Several  errors  could  contribute  to  the  results  at  the  locations  v^ere  the  GCI  did  not  lie  in 
the  asyiiq>totic  range.  First,  Roache  (1993)  mentioned  uang  geostatistical  methods  to  generate 
grid-block  property  variations  with  specified  statistical  parameters,  and  the  question  of  Tvbether 
the  propoty  fields  should  change  as  a  result  of  a  changing  grid  structure,  or  be  fixed.  These  head 
solutions  assumed  changes  in  the  geological  variables  as  the  grid  became  coarser.  The  hydrauEc 
conductivity,  for  exarxple,  was  re-kriged  for  every  changing  grid  size.  This  mediod  was  much 
eaaer  to  simulate,  but  probabty  contributed  to  the  different  head  sohiticnis  bemg  out  of  the 
atyu^totic  range.  Secondty,  the  numerical  convergence  criterion  in  the  flow  model  was  constant 
for  all  three  solutions.  This  may  cause  over-confidence  in  tire  coarser  grid  solution.  Lastty, 
machine  precision  may  have  caused  errors  to  propagate  throughout  the  emulations.  Future 
sirtrulations  ^ould  use  double  precision  in  the  modeling  to  reduce  machr  ‘^related  errors. 
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Since  none  of  the  solutions  were  in  the  asynq)totic  range,  an  error  based  on  the  grid 
pacing  calculations  could  not  be  coni^uted.  The  calculated  GCIs  were  virtually  meaningless. 
Future  emulations  should  use  fixed  geological  property  fields,  eq>ecially  hydiauHc  conductivity. 
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SECTION  IX 
CONCLUSION 


A  numerical  model  of  the  MADE-2  e7q)ernneat  was  created  to  amulate  groundwater  flow 
and  contaminant  transport  using  the  pubfic  domain  programs  MODFLOW  and  MT3D.  The 
following  conclusions  were  reached. 

1.  The  use  of  public  domain  software  on  a  Pentiomrclass  personal  conoputer  to  model  an 

experiment  such  as  MADE-2  is  very  practical  MODFLOW,  \^Mch  solved  the  groundwater  flow 

equation  for  the  hydraulic  heads,  did  a  good  job  in  matching  tiie  observed  phenomena.  The 

maximum  average  RMS  difference  for  any  stmnlation  was  0.23  meters.  This  translates  to  an 

•  •  2 

average  23  cm  difiference  ovmnll  b^weoi  heads  of  the  observed  and  .simulated  within  a  33000  m‘ 
domain. 

The  contaminant  tran^ort  model,  solved  by  MT3D,  was  less  successfiiL  The  final 
agreement  of  the  model  with  the  observations  was  good,  but  the  necessary  assumptions 
concerning  horizontal  anisotropy  have  not  been  justified.  Of  course  tire  simulated  plume  of  any 
contaminant  does  not  match  exactfy  with  the  observed.  Small  scale  heterogeneities  in  the  field  are 
impoRsihle  to  reproduce  using  geostatistical  techniques  with  a  finite  data  set.  Additional  hydraulic 
conductivity  measurements,  however,  would  probabfy  not  change  the  resuhs  dramatically. 

2.  To  accurately  modd  the  tran^ort  of  trithimj  a  horizontal  anisotropy  foctor  idiich 
varied  by  postion  had  to  be  introduced.  F^arHer  modeling  efibrts,  which  assumed  isotropic 
conditions,  did  not  allow  the  phune  to  migrate  as  for  as  the  observed  phime.  Punqi  test  AT-2 
suggested  horizontal  anisotropy  with  the  major  princ^al  axis  35**  west  of  north,  in  Boggs  et  aL 
(1992),  Young  (1995),  and  other  pikers,  a  meandering  abandoned  stream  channel  was 
hypothesized,  whidi  cut  throu^  the  MADE-2  site  at  qiproxiniatefy  30°  east  of  north.  Realistic 
phune  predictions  were  achieved  on]^  by  agsiiming  the  chaimel  did  in  ftict  exist  and  was  located  as 
predicted.  However,  the  assumed  princqial  axes  were  not  aligned  with  the  diaimel  axis  or  with 


164 


the  AT>2  drawdown  ell^ses.  The  assumed  major  pimcq)al  axis  was  aligned  with  the  Y  axis,  12*’ 
west  of  north.  In  addition,  the  assumed  anisotropy  &ctor  inq>£ed  higher  conductivities  than  were 
measured  with  the  borehole  flow  meter,  and  was  hi^er  than  any  found  in  the  literature.  In  short, 
the  horizontal  anisotropy  needed  to  match  the  observed  plumes  is  not  fiiOy  supported  by  the  field 
data. 


3.  Knowing  the  plumes,  a  priori,  was  essential  in  accurately  modeling  the  experiment. 
Without  such  knm\iedge  it  would  have  been  inq>ossible  to  estimate  the  values  of  diq>ersion  or 
anisotropy  needed  to  push  the  plume’s  tail  into  the  far  field.  This  is  very  discouraging,  as  the 
purpose  was  to  predict  the  plume’s  movement,  and  not  to  just  rqtroduce  it.  Future  experiments 
at  the  MADE  site  will  at  least  know  something  about  the  values  of  di^^on  or  anisotropy. 
However,  modelers  of  other  extreme^  heterogeneous  ates  will  not  be  so  fortunate.  The  tools 
used  in  this  study  do  not  provide  a  true  predictive  cq>abiiity. 
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APPENDIX  A 


KRICSED  OBSERVED  HEAD  SURVEYS 
(See  Figure  13  for  June  19,1990,  and  Hgure  14  for  March  18, 1991.) 
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Upper  Screened  Heads  of  October 
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Figure 


5,  1990 


upper  Screened  Heads  of  January 
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X(m) 


Upper  Screened  Heads  of  August 
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APPENDIX  B 


KRIGED  TRANSMISSIVITY  AND  LEAKANCE  FOR  M2-8-8 

Based  on  bordiole  flowmeter  data  throng  K81 
(See  Figure  19  for  Layer  3.) 
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[m^  /d]  of  Layer  4. 
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FiguieB-15.  Leakance  [1/d]  of  Layer  8. 
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APPENDIX  C 

CHANGES  TO  MODFLOW  SOURCE  CODE  FOR  MF95 


These  are  the  changes  to  the  MODFLOW  source  code  to  produce  MF95.FOR; 


1)  lu  the  BCF2AL  subroutine,  space  allocation  needs  to  be  larger. 


Qiange: 

LCTRPY  =  ISUM 

ISUM  =  ISUM+NLAY 

To: 

LCTRPY  =  ISUM 
ISUM  =  ISUM+ISIZ 

2)  In  the  BCF2RP  subroutine. 

Change: 

DIMENSION  HNEW(NODES),  SCl(NODES),  HY(NODES),  CR(NODES), 
1  CC(NODES),  CV(NODES),  ANAME(6,11),  DELR(NCOL), 

1  DELC(NROW),  BOT(NODESX  TOP(NODES),  SC2(NODES), 

1  1RPY(NLAY),  IBOUND(NODES),  WETDRY(NODES), 

1  CVWD(NODES) 

To: 

DIMENSION  HNEW(NODES),  SCl(NODES),  HY(NODES),  CR(NODES), 
1  CC(NODES),  CV(NODES),  ANAME(6,11),  DELR(NCOL), 

1  DELC(NROW),  BOT(NODES),  T0P(N0DES),  SC2(N0DESX 

1  TRPY(NODES),  IBOUND(NODES),  WETDRY(NODES), 

1  CVWD(NODES) 


and  change: 

CALL  U1DREL(1RPY,ANAME(1,8),NLAY,IN,I0UT) 
CALL  UlDRELpELR,ANAME(l,9),NCOL,IN,IOUT) 
CALL  U1DRELPELC,ANAME(1, 10),NROW,IN,IOlfr) 
C 

C2 - READ  ALL  PARAMETERS  FOR  EACH  LAYER 

KT=0 

KB=0 
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DO  200  K=1,NLAY 
KK=K 
C 

C2A— FIND  ADDRESS  OF  EACH  LAYER  IN  THREE  DBVIENSION  ARRAYS. 
IF(LAYCON(K).EQ.l  .OR  LAYCON(K).EQ.3)  KB=KB+1 
IF(LAYCON(K).EQ.2  .OR  LAYCON(K).EQ.3)  KT=KT+1 
LOC=l-KK-l)*ND 
LOCB=l+(KB-l)*ND 
LOCT=l-KKT-l)*ND 
C 

C2B - READ  PRIMARY  STORAGE  COEEFICIENT  INTO  ARRAY  SCI  IF  TRANSIENT 

IF(ISS.EQ.O) 

CALL  U2DREL(SCl(LOC)^AME(  1,  l),NROW,NCOL^IN,IOUT) 


to: 

CALL  UlDREL(DELR^AME(l,9),NCOL,IN,IOirr) 

CALL  U1DREL(DELC,ANAME(1,10),NROW,IN,IOUT) 

C 

C2 - READ  ALL  PARAMETERS  FOR  EACH  LAYER 

KT=0 

KB=0 

DO200K=l,NLAY 

KK=K 

C 

C2A - FIND  ADDRESS  OF  EACH  LAYER  IN  THREE  DIMENSION  ARRAYS. 

IF(LAYCON(K).EQ.l  OR  LAYCON(K).EQ.3)  KB=KB+1 
IF(LAYCON(K).EQ.2  .OR  LAYCON(K).EQ.3)  KT=KT+1 
LOC=l-KK-l)*ND 
LOCB=l-KKB-l)*ND 
LOCT=l-KKT-l)*ND 
C 

C2B - READ  PRIMARY  STORACffi  COEFFICIENT  INTO  ARRAY  SCI  IF  TRANSIENT 

CALL  U2DREL(lRPY(LOC)^AME(l,8)^OW,NCOL^IN,IOUT) 
ffaSSLQ.O) 

CALL  U2DREL(SCl(LOC),ANAME(l,l),NROW,NCOL,KK4N,IOUT) 
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*Note:  The  new  TRPY  vaiiable  is  to  be  read  by  the  U2DREL  array  reader  before  eadh  SFl  array 
(if  transient,  or  transmissivity  if  steady*state),  as  opposed  to  being  read  by  the  UIDRH-  array 
reader  before  DELR. 


3)  ha  BCF2FM,  SBCFIC,  SBCF2H,  and  SBCF2N,  changes  need  to  be  made  by  changing  the 
dimensions  of  TRPY  (NLAY)  to  TRPY  (NCOL,  NROW,NLAY). 


4)  The  last  changes  are  to  be  made  to  the  SBCFIC  subroutine. 
Change: 


YX=TRPY(K) 

C 

Cl - ^FOR  EACH  CELL  CALCULATE  BRANCH  CONDUCTANCES  FROM  THAT  CELL 

Cl - TO  THE  ONE  ON  THE  RICaiT  AND  THE  ONE  IN  FRONT. 

DO40I=l,NROW 

DO40J=l,NCOL 

T1=CC(J,I3C) 

to: 


C  YX=TRPY(K) 

C 

Cl - FOR  EACH  CELL  CALCULATE  BRANCH  CONDUCTANCES  FROM  THAT  CELL 

Cl - TO  THE  ONE  ON  THE  RIGHT  AND  THE  ONE  IN  FRONT. 

DO40I=l,NROW 

DO40J=l,NCOL 

T1=CC(J4,K) 

YX  =  TRPY(TO 


This  is  the  aid  of  the  changes  to  MODFLOW. 
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